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ABSTRACT 
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^c^^. We describe calculations for the formation of icy planets and debris disks 

at 30-150 AU around 1-3 Mq stars. Debris disk formation coincides with the 
'^ . formation of planetary systems. As protoplanets grow, they stir leftover plan- 

"^ . etesimals to larere velocities. A cascade of collisions then errinds the leftovers 

to dust, forming an observable debris disk. Stellar lifetimes and the collisional 
cascade limit the growth of protoplanets. The maximum radius of icy planets, 
f'max ~ 1750 km, is remarkably independent of initial disk mass, stellar mass, and 



O, 

00 , stellar age. These objects contain < 3%-4% of the initial mass in solid material. 

Collisional cascades produce debris disks with maximum luminosity ~ 2 x 10^^ 

times the stellar luminosity. The peak 24 yum excess varies from ~ 1% times the 

stellar photospheric flux for 1 M© stars to ~ 50 times the stellar photospheric 

flux for 3 Mq stars. The peak 70-850 /xm excesses are ~ 30-100 times the stellar 

photospheric flux. For all stars, the 24-160 /xm excesses rise at stellar ages of 

5-20 Myr, peak at 10-50 Myr, and then decline. The decline is roughly a power 

law, / oc t~" with n ^ 0.6-1.0. This predicted evolution agrees with published 

observations of A-type and solar-type stars. The observed far-IR color evolution 

of A-type stars also matches model predictions. 

Subject headings: planetary systems - solar system: formation - stars: formation 
- circumstellar matter - infrared: stars 
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INTRODUCTION 



During the past 25 years, observations from IRAS, ISO, and Spitzer have revealed sub- 
stantial m id-infrared (mid-IR) excesses associated with hundreds of normal main sequence 



stars (e.g..lBackman fc 



2006; Moor et al. 



2006 



'aresce 



1993 



Habing et al. II2OOII : iRieke et al. II2005I : iBrvden et al 



Rhee et al. I l2007al ) . Curren t samples include s tars with spectral 



types A-M and ages ~ 5 Myr to ~ 10 Gyr (e.g., 



Rieke et al. 2005 



Beichman et al. 2006; Su et al. 



Chen et al. 



2006 



2005; 



Kim et al. 2005 



Hillenbrand et al. II2OO8I ). AL 



though binary stars and single stars in dense clusters and in the field are roughly equally likely 



2006; Currie et al. 



to have IR excesses (Stauffer et al. 112005: 



2007a 



Siegler et al. 



Su et al. ll2006l:lBrvden et al. I bood: ICorlova et aL 



emission declines from ~ 30%-40% for A- type stars (ISu et al 



solar-type stars (JGreaves &: Wyatt I l2003l ; iTrilling et al. 



200?! : iTrilling et aUboojVtli e frequency of 



2008 



excess 



2006h to ~ 10%- 20% for 
Mever et al. Iboosh . Thus, 



this phenomenon is common among main sequence stars and may depend on stellar mass. 
Hi gh quality images denaonstrate that dust orbiting the central star produces the ex- 



cesses (ISmith fc Terrile Ill984l : iBrandeker et al. 112004 : IStapelfeldt et al. 112004 : iKalas 112005 



Meyer et al. 1 120071 ) . In j3 Pic and AU Mic, t he dust is in a geometrically thin, edge-on 



disk with an outer radius of a ~ 200-1000 AU (ISmith fc Terrile 



Golimowski et al. 1993; Kalas et al. 2004; Liu 



2004 



1984 



Augereau fc Beust 



Telesco et al. 111988 
2006h . In these 



disks, the small scal e height of the dust, H/a ~ 0.1, is consistent with material in roughly 



circular orbits (e.g.. iBackman fc Paresce Ill993l : iKalas et al. II2004I . an d references therein) . 



Although broad tori of dust are visible in many other systems (e.g. iGreaves et al 



1998 



Augereau et al. Ill999l : iHoUand et al. II2003I : iKalas et al. II2OO6I : ISu et al. II2OO6I: Fitzgerald et al. 



2007 ) , narrow rings of dust produce the emission in a PsA and HR 4796A (IJavawardhana et al 



Schneider et al. Ill999l : ICreaves et al. II2OOOI : iTelesco et al. II2OOOI : iKalas et al. Il2005[ ). 



1998 



For systems with face-on rings and tori, the total emission constrains the scale height, H/a ~ 
0.1. Thus, the dust in these systems is as highly flattened as the structures in /? Pic and AU 
Mic. 

Broadband spectral en ergy distributions constrain the luminosity, size, t emperature. 



and total mass of the dust (IBackman fc Paresce I Il993l : iLagrange et al. I I2OOOI : iDent et al 



2OOOI : IWolf &: Hillenbrand I l2003l ) . Some stars have exces s es from grains plausibly associ- 



2007bl : iRhee et al. Il2007bl : iMever et al. 



ated with the terrestrial z o ne (e.g.. iBeichman et al. II2005I: lAbsil et al. I l2006l : I Currie et al 



2008 



Lisse et al. I l2007d ) . Optical ly thin emission 



from cooler materia l with temperature T ~ 2 0-150 K is more typical (e.g. ISu et al. 112006 



Trilling et al. I l2008l : [Hillenbrand et al. II2008I ). For systems with submm observ ations, the 
measured fluxes suggest grains with sizes ~ 1 um-1 cm and total mass ~ .01 M^ (ILiu et al. 



2OO4J : iNaiita fc WiUiams I l2005l : IChen et al. I I2OO6I : IWilhams fc Andrews 1 120061 ) . The grains 
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have compositions s imilar to dust in the asteroid belt, comets, or the trans- Neptu nian region 
of the Solar System (JGriin et al. Ill995l : iBrownlee et al. 1119971 : iLisse et al. ll2007al Jbl) . Because 
the dust mass in these systems lies between the initial mass of solids in protostellar disks 
100-1000 M®; iNatta et al.lboool: [Andrews fc WiUiams Iboosi) and the dust mass in th e 



Solar System (< 10"^ M®; Iflahn et al. II2OO2I : llanderaf etliLlbood : iNesvornv et al. Ibooeh 



the dusty structures in the se systems are often called 'debris disks' (JBackman fc Paresce 
I993I : iLagrange et al. ll200o[ ). 



In addition to the dust properties, several other observations suggest plausible links be- 
tween debris disks and the formation of plai ietary systems. Obse rvations of A-type stars sug- 
gest a 'rise and fall' of debris disk emission (jCurrie et al. Il2008al ). with a clear increase in the 
typical 24 /im excess at 5-10 Myr, a peak at 10-20 Myr, and a decline for t > 20-30 Myr. The 
rise in debris disk emission roughl y coincides with the disappea. r ance of optically thick emis 



Hernandez et al. 



sion from protostellar disks fe.g. iHaisch. Lada. fc Lada 1 1200 It ISicilia-Aguilar et al. 1 12006 



20061: ICurrie et al. ll2007al Jd'). The broad plateau occurs at a time when ra - 



dio metric dating (JYin et al. II2002I ) and theory (jChambers II2OOII : iKenyon fc Bromley II2006I ) 



suggest the Earth contained ~ 90% of its final mass. The decline of dusty debris around A- 
type stars is at ro ughly the same time as a gradual decrease in the cratering rate of objects in 



the Solar Syste m fISwindle lll993l : lMelosh. Vickerv. fc Tonks lll993l : IWadhwa fc Russell II2OOO 



Koeberl II2003I ). These results suggest that the evolution of dust in debris disks parallels the 
evolution of larger solid objects in the Solar System. 

Simple physical arguments also link debris disks with the formation of planetary systems. 
Because radiation removes 1-100 /im grains on timescales shorter than the stellar age, some 
process replenishes the dust. To maintain the observed dust masses for long timescales, 
normal stars must have a large reservoir, ~ 10-100 M®, of unseen objects that continuously 
collide at large velocities and fragment into smaller objects. Remnant material from planet 
formation satisfies both needs. The growth of 1000 km or larger planets in a disk of small 
grains naturally leaves behind an e nsemble of 'leftover' 1-10 km 'planetesimals' on eccentric 
orbits (JKenyon &: Bromley I l2004al 2004b). For a mass of ~ 10 -100 M^ in leftovers, high 



velocity collisions produce enough dust for most debris disks (e.g.. lBackman fc Paresce Ill993 



Habing et al. II2OOII : iKenyon fc Bromley Il2004bl ). If this interpretation is correct, debris disks 
provide conclusive evidence for the formation of Pluto-mass or larger planets around many, 
if not most, main sequence stars. 

In addition to these considerations, numerical calculations suggest that an evolving 
swarm of 1-10 km planetesimals explains several observed trends in the properties of de 



bris dis ks. Starting with an ensemble of < 1 km-sized planetesimals, IKenyon &: Bromley 



(j2002bl . 2004a, b, 2005) show that collisions and mergers form 500-1000 km-sized objects 
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in 1-50 Myr. These protoplanets stir up leftover planetesimals along their orbits. De- 
structive collisions among the leftovers then produce a collisional cascade - where col- 
lisions gradua lly grind large objects into smaller ones - alorig with copious amounts o f 
dust (see also Wi^ 



Dominik fc Decin 



lianas fc Wetherill 



1994 



Durda fc Dermott 



1997 



Quillenet al. 1120071 ). 



J2003h . IWvatt et al.1 J2007ah . and llohne et al.1 (|2008|) show that colli 



sional evolution in a belt of high velocity planetesi mals naturally produces a d ust luminosity 
that declines roughly inversely with time (see also lKenyon fc Bromley Il2002al 2004b, 2005), 
explaining the observed time evolutio n - L^ (x t~", with n ~ P-5~l - suggested by re- 
cent observations of A-type stars (e.g.. iKalas Ill998l: iHabing et al. Il200ll : iDecin et al. 112003 



Greaves fc Wvatt II2003I : iRieke et al. II2005I : iRhee et al. 



2007af) . To ac c ount fo r the large ob 



served range of IR excesses among stars of similar ages, IWyatt et al. I (l2007al ) propose belts 
with a range of initi al masses and semimajor ax es, as suggested from submm observations 
of protostellar disks ([Andrews &: Williams II2005I . 2007b). 



Here, we continue to explore the evolution of dusty debris arising from planet formation 
in a disk of icy planetesimals. Our suite of calculations for disks at 30-150 AU around 1- 
3 M0 stars yields robust predictions for the maximum sizes of icy planets as a function of 
semimajor axis and stellar age. Results for the long-term evolution of IR excesses account 
for many fundamental aspects of the data. T hese calculations are t he first to explain the 'rise 
and fall of debris disks' around A-type stars (ICurrie et al. Il2008al ) and the apparent peak i n 
the 70-160 /im excesses of G-type stars with ages of ~ 100 Myr (JHillenbrand et al. II2008I ). 



Comparisons between our models and current observations suggest that the minimum stable 
grain size and the slope of the IR emissivity law are critical parameters. 

The models suggest a set of further critical observations. Spatially resolved images 
of debris disks around A-type and solar-type stars can improve our understanding of the 
minimum stable grain size. Larger samples of debris disks with high quality submm data 
from ALMA, Herschel, and SOFIA can place better constraints on the slope of the emissivity 
law. Together, these data can test our predictions for the time evolution of debris disk 
emission around 1-3 M© stars and provide input for more complete calculations that include 
the formation and dynamical evolution of giant planets. 

We outline our model in §2. We describe results for the formation of icy planets in §3 and 
the evolution of debris disks in §4. After discussing several applications of our calculations 
in §5, we conclude with a brief summary in §6. 
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THE MODEL 



Kenvon fc Bromlev I (1200 ll . 2002a, 2004a, 2004c) and iBromlev fc Kenvon I (120061) de- 



scribe our hy brid multiannulus nu r aerica l model for planetesimal growth. iKenvon fc Luu 
Jl998l . 1999). iKenvon fc Bromlevl J200ll . 2002a), and iBromlev fc Kenvon l J2006h compare 



results with analytical and numerical calculations. We adopt the ISafronov I (Il969l ) statistical 
approach to calculate the co llisional evolution of an ensemble of planetesimals i n orbit around 



a star of mass TkL (see also ISpaute et al. 



20061 : iThebault fc Augereau 



2007 



1991 



Weidenschilling et al. 1997: Krivov et al. 



Lohne et al. I l2008l ) . The model grid contains A^ con- 



centric annuli with widths 5ai centered at semimajor axes Oj. Calculations begin with a 
differential mass distribution n{mik) of objects with horizontal and vertical velocities /ijfc(t) 
and Vik{t) relative to a circular orbit. The horizontal velocity is related to the orbital ec- 
centricity, e'^i^it) = 1.6 {hik{t)/VK,iY, where Vx^i is the circular orbital velocity in annulus i. 
The orbital inclination depends on the vertical velocity, ijj^it) = sm~^ {2{vik(t) /VK,iY) ■ 

The mass and velocity distributions evolve in time due to inelastic collisions, drag forces, 
and gravitational forces. For inelastic collisions, we solve the coagulation equations for a 
particle in mass batch k of annulus i colliding with another particle in mass batch I of 
annulus j, 



6n. 



i'k' 



St [^ijkiAijkiTiiknji — ni'k'Ai/jk'iriji] + Sriiiyj — Suiiy^gd (1) 

SMi'k' = St mi'k' [^ijkiAijkinikTiji — ni/k'Ai/jk'inji] + 6Mi'k'j — SMi/k'^gd (2) 

where t is time, Mj/fe/ is the total mass in mass bin k' in annulus i', Aijki is the cross-section, 
^ijki = 1/2 for i = j and k = I, and eijki = 1 for k ^ I and any z,j. The terms in these 
equations represent (i) mergers oirriik and rriji into a body of mass rrii'k' = rnik + rnji — rrie^ijki, 
(ii) loss of mj'fc' through mergers with other objects, (iii) ad dition of mass from debris 
prod uced by the collisions o f other objects (JKenyon fc Luu Ill999l ). and (iv) loss of mass by gas 
drag (JKenyon fc Luu Ill998l ). In each equation, the second term includes the possibility that a 



collis ion can produce debris but no merger (rebounds; see lDavis et al. lll985l : lKenyon &: Luu 
19991 . and references therein). 



The collision cross-section is 

A-ijkl = Oicoll 



Vijki Fg^ijki {rik + Tji) 



(3) 



^4 Hijki {aij) {Aaij)^ 

where acou is a constant (IWetherill fc Stewart lll993l : lKenyon fc Luu Ill998l ). H^ki = J 2 {vf^. + v'ji)/ {Vlij) 
is the mutual scale height, {aij) and (Aaij) are the average heliocentric distance and width for 
the two interacting annuli, (Omegaij) is the average angular velocity, Vijki is the relative par- 
ticle velocity, Fgijki is the gravitational focusing factor, and r^k and Tji are the particle radii. 
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We adopt the piecewise analytic approximation of ISpaute et al. I (Il99ll ) fo r the gravitational 



focusi ng factor in the dispersion regime and the c oUisional cross-sections of iGreenberg et al 



( ll99ll ) in the shear-dominated regi me (see also iGreenzweig &: Lissauer 



more details of this algorit hm, see iKenvon &: Luu 



Kenyon fc Bromley I (J2004al ) , and iBromley fc Kenyon I (120061 ) . 



Jl998h 



1990 



199 2). Fo r 



Kenyon &: Bromley I (l2002al ). 



To choose among possible collision outcomes, we use an energy-scaling algorithm. If Q'^ 
is the collision energy needed to eject half the mass of a pair of colliding planetesimals and 
Qc is the center of mass collision energy, the mass of the ejecta is 



m. 



e,ijkl 



0.5 {rriik + niji) i — ^ j 



(4) 



where rrtik and ttiji are the masses of the colliding planetesimals. This approach allows us to 
derive ejected ma sses for catastrophic collisions with Qr ^ Q^ and cratering collisions wit h 



Qc < Q*d (see also IWetherill fc Stewart Ill993l : IStern fc Colwell 1 19971: iKenvon fc Luu Ill999h . 



Consistent with N-body simulations of collision outcomes (e.g. 
set 



Benz &: Asphaug Ill999l ). we 



A 



Qbr^ji^l + QgPgl^ijkl 



(5) 



where rijki is the radius of a merged object with mass rriik + rriji, pg is the mass density of a 
planetesimal, Qhr^*' is the bulk component of the binding energy, and QgPgV^^ is the gravity 
component of the binding energy. 



Kenyon fc Bromley I (120051 ) and iKenyon et al. I (120081 ) describe how collisional evolution 
depends on various choices for Q^. For icy objects, detailed numerical collision simulations 



< 



10^ erg cm ^'' 



K. 



-0-5 < 0. < 0, 



A^ 



yield Qb 

g~^, and Pg ~ 1-2 (e.g., iBenz fc Asphaug Ill999l : iLeinhardt et al. 



1-2 g cm ^, Qg < 1-2 erg cm^ ^^ 
20081 ). Calculations for 



the breakup of comet Sho emaker-Levy 9 suggest a smaller component of the bulk strength 



QbT^'' ~ 10 erg g (e.g.. lAsphaug fc Benz Ill996l ). which yields smaller disruption energies 



for smaller objects. Because ne arly all models for collisional disruption yield sim ilar results 



for objects with r > 1 km (e.g., iKenyon &: Bromley Il2004d : iKenyon et al. II2008I ). collisional 



evolution is relatively independent of these uncertainties as planetesimals grow into larger 
objects. Thus, we choose standard values - Qg = 1.5 erg cm^-^^ g~^, pg = 1.5 g cm~^, and 
f3g = 1.25 - for the gravity component of Q"^. To check how the evolution of the small 
planetesimals depends on Q^, we conside r a broad range in the bulk component of t he 



2005 



Kenvon et al. 1120081 ) 



strength, Qb = 1^10^ e^K K ^ with /?fe = (jPan fc Sari I l2005l : IKenvon fc Bromlev 112004c 



To compute velocity evolution, we include collisional damping from inelastic collisions, 
gas drag, and gravitational interactions. Our equations for the evolution of the velocity 
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dispersion are 



dhl 



"'"'iriAk 



dt dt 

for the horizontal component and 






dt 



dt 



dh"^ 

'-''"'sr,ik 

dt 



dt 



dt 






dv^ 
dt 



(6) 



(7) 



for the vertical component, where the subscripts refer to the contributions from collisional 
damping ('in'), gas drag ('gd'), and long-range ('Ir') and short-range ('sr') gravitational 
interactions. 



For collisional damping, we adopt 



and 



dhf^,ik 
dt 



d'^inAk 



3 -^* ^ — ''•max 

^ J2 ^ ^"^ji^^i ~ ^^khj 

j=0 1=0 






ik 



{niik + mji)hl^) leiPijki) 



j=Nl=l„ 



tvii 



dt 






P! 



[rrijiv'^i - rriikVi^ - (rriik + 'mji)Vii.) Ii{Pijki) 



(9) 
where Cin = acoii fijki eiju Pji Vijki Fg^ijM (r^fc + rjif, Pf-^i = {tf^ + i]M {e% + e|, ), and 



Pu is the volume density of planetesimals with mass rriji in annulus j (lOhtsuki I Il992 



Wetherill fc Stewart 



1993). In the second summation, L, 



k when i = j; Ir. 



M 



when i ^ i (see also iKenyon fc Luu Ill998l . 1999). We add a term, fjjku to treat the overlap 
between adjacent zones; fijki = 1 when i = j and fijki < 1 when i j^ } (jKenyon &: Bromley 



20011). The integrals I^ and Ij are elliptic integrals described in previous p ublications 
JWetherill fc Stewart 1^9931 : Istewart fc Idalboool : lOhtsuki. Stewart, fc Ida]l2002h . 



For velocity damping from gas drag, we follow IWetherill &: Stewart I (119931 ) and write 

(^'T'gd,ik 



dt 



Pik ^ Pgas^ gas^ik) 



2mik 



and 



dv, 



gd,ik 



dt 



-(1-Afc) T^PgasVgasr. 



2m. 



ik 



gas ik^ 



(10) 



(11) 



where Cd = 0.5 is the drag coefficient, [3ik = hik/{hi^ + w^) ' , Pgas is the gas density, r] is 



the relative gas velocity, and Vgn.!, = (VikJVik + v) 



see 



Adachi et al. Ill976l : IWeidenschilling Ill977b 



^/^ is the mean relative vel ocity of the gas 
Wetherill fc Stewart I[l993l l 



For gravitational interactions, we compute long-range stirri ng from distant oligarchs 
( Weidenschilling Ill989l ) and short-range stirring from the swarm ( Ohtsuki. Stewart, fc Ida 



20021 ). The long-range stirring only has a horizontal component, 

G PjiMji /tan" {Hijki/Dmin) tan" (Hij^i/Dmax) 



"'"'IrAk 



dt 



j=N l=M 

i=l 1=1 



(%) 



D, 



D„ 



(12) 



for continuum objects and 



dhj^ik 
dt 



j=N l=M 

EE 



G^ 



^lr,e'^jl 



^ ^^ TiVta \ i5a? + 0.5if 
I 1=1 \^ i 



2\2 



(13) 



for individual oligarchs, where Xijki is the fraction of objects with mass rriik in annulus i 
that approach no closer than 2.4 Rh of the objects with mass mi in annulus j, -D^j^ = 



max(2.4i?H, 1.6(/i?;. + /^^^)^^^), -Dmo,.T: = 0-5 ma xfw,:^, Wi;), (^a 



la,; — a,- 1, C] 



^ire — 5-9 (sGe also lKenyon fc Bromley 1120011 ). 



'J I ) ^lr,e 



23.5, and 



For short-range gravitational interactions, the stirring depends on the ratio of the rela- 
tive collision velocity to the mutual Hill velocity. 



Vh 



(%) {aij) [{mik + m,i)/3M,Y/^ . 



(14) 



In the high velocity regime, the collision velocity exceeds the Hill velocity. Statistical solu 
tions to the Fokker-P l anck equation then yield accurate estimates for the stirring rates (e.g., 



Hornung et al. Ill985l : IWetherill fc Stewart Ill993l : IStewart k Ida II2OOOI : iKenvon fc Bromlev 



2OOII) . At low velocities, n-body calculations provide good estimates. We follow lOhtsuki. Stewart, fc Ida 
(J2OO2I ) and write the stirring as the sum of rates in the two regimes: 



and 



^"'sr,ik 

dt 



dv"^ 

"-^srAk 



"'"'high,ik 
Jt 



"'^high,ik 



+ 



dh"^ 

'-^"'low,ik 
Jt 



dt 



dt 



di? 

"'^low,ik 

dt 



(15) 



(16) 



wher e the subscripts 'high' and 'low' indicate the velocity regime (e.g., Eq. (25) of lOhtsuki. Stewart, fc Ida 
2002h . 



In the high velocity regime, the stirring is (e.g. lStewart fc Ida ll2000l : lKenyon fc Bromley 



2OOID 



^^2 j=Nl=M 

high,tk = ^ ^ f^jf^i Ghigh {{hl^ + h]i) rriji Je{(3ijki) + IA {rrijih^i-mikhl,}) He{f3ijki)) (17) 



dt 



j=i 1=1 



-9- 



and 

J 2 j=Nl=M 



dt 



j=l 1=1 



^high 
f^iik 



^ X] f'^''^ ~iW^ ii^ik+t^'ji) ^ji Jz{(iijki) + IA {mjivl-m,kVik) H^{(3ijki)) , (18) 



ijkl 



where fijki the fraction of objects with mass rriik in annulus i that approach within 2.4 Rh 
of the objects with mass rriji in annulus j and Chigh = 0.28 A\ G^ Pii/((^L + ^j/)'^^^)- ^^ ^^e 
expression for Chigh, ^a = In (A^ + 1) and 



A 



'0.19 (hl + hl + 1.25 ivl + v],)) ivl + v],y/'^ 



(19) 



-'H 



The functions He, H^, Je, and Jz are definite integrals defined in Stewart & Ida (2000). 



In the low velocity regime, the evolution equations are (jOhtsuki. Stewart, fc Idall2002l ): 

(20) 



iow,ik = ^ ^ f..f^i Cion, {rriji xi + {mjih'^ji - rriikhli?) Xs) 



and 



dt 



di? 



j=l 1=1 

j=Nl=M 



dt 



X^ X^ fijki Clow (rriji X2 + ijrijih^i - rriikhlj.) Xi) 
j=i 1=1 



where the x's are simple functions of the Hill radius 

rH,ijH = a[{2 {mik + mji))/3M^Y^^ . 



(21) 



(22) 



and the normalized eccentricity and inclination (jOhtsuki. Stewart. &: Ida 1120021 . see also Ida 
1990; Ida & Makino 1992). For the low velocity limit of the horizontal velocity 



Xi - 73 Ci rfjij^i 



and 



X2 = C2{Atl + 0.2Bl{eltiy/')rj,^,^,i. 
For the low velocity limit of the vertical velocity 

X3 = 10 C, ~e% r%^,,Jihl + hi) 

and 

X4 = 10 Cs i% rjj,ijki/{h'^ik + h^) ■ 



(23) 
(24) 

(25) 
(26) 



Here, Qow = 0.625 (atj) pji Hijki VlJ{mik + m^if , l\ = (e^^ + e]i)/r]j^-^i, and i% = 
{i% + i'j))/rffij^.i- The constan ts Ci, C2, and C3 are identical to those in Eq. (26) of 
Ohtsuki. Stewart, fc Ida! J2002h . 
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Several tests indicate that these expressions provide an accurate trea tment of velocity 



evolut ion for planetesimals in the high and low velocity regimes. Figs. 5 -7 of lOhtsuki. Stewart, fc Ida 
( 20021) show comparisons with results from n- body simulations (see al so Ida lllQQOl : Ida fc Makino 



9921) . Our simulations confirm this analysis. IWeidenschilling et al. I ( 19971) and Kenvon fc Bromley 



(2001) compare variants of this formalism with other n-body calculations. iGoldreich. Lithwick. fc Sari 
(J2004j ) demonstrate that our numerical results agree with analytic estimates. 



To follow the evolution of the most massive objects more accurately, our code includes 
an n-body algorithm. When objects have masses exceeding a 'promotion mass' rripro, we 
promote them into an n-body code that directly solves their orbits. The n-body code incor- 
porates algorithms to allow collisions among n-bo dies and interactions between n-bodies and 
coagulation particles. iBromley fc Kenyon I (120061 ) describe this code in detail. Because dy- 
namical interactions among large oligarchs are rare and occur at late stages in the evolution, 
we set m„rn = 10^^ g. We describe several test calculations with smaller mpm in ^3.3. 



"pro 



"pro 



To treat the time evolution of the gas volume density pgas, we use a simple nebular model 
with gas surface density T,gasia,t) = T,gas,oa~^^'^e~*^*^''% gas-to-sohds ratio T.gas,oi<^)/^io.) = 



100 - where S is the surface density of sohds, and scale height Hgas{a) = Hgasflia/aoY-^'^^ 



gasy-") ^^gas,\ 

(JKenyon fc Hartmann 1 119871 ). To approximate gas removal on a timescale tgas-, the gas 
density declines exponentially with time. We set tgas = 10 Myr. During the early stages 
of calculations at 30-150 AU, velocity damping is important for particles with r < 100 m. 
However, particle losses from gas drag are small, ~ 1% or less of the initial mass. By the time 
viscous stirring dominates the velocity evolution of small object s, the gas disk has dispersed . 
Inward drift and velocity damping are then negligible (see also IWetherill fc Stewart Ill993l ). 



The initial c onditions for these calculation s are appropriate for a disk with an age of 



< 1-2 Myr (e .g. IPuUemond fc Dominik I l2005l : iNomura fc Nakagawa I l2006l : ICiesla 1 120071 : 



Garaud 1120071 ). We consider systems of A^ annuli in disks with Oj = 30-150 AU and Sai/ai 
= 0.025. The disk is composed of small planetesimals with radii of ~ 1-1000 m and an 
initial mass distribution rtiimik) in each annulus. The mass ratio between adjacent bins is 
5 = rriik^i/mik = 1.4-2. At the start of the calculations, each bin has the same total mass, 
eccentricity cq = 1 — 3 x 10~^, and inclination iq = eo/2. We assume a power law variation 
of the initial surface density of solid material with semimajor axis. 



So(M. 



(27) 



where Xm is a scaling factor. For a 1 Mq central star, models with Sq ~ 0.1-0.2 g cm ^ 
at Oq = 30 AU h ave a mass in icy solids comparable to the minimum mass solar nebula 



(MMSN hereafter; IWeidenschilling Ill977al : iHavashi 



disks surrounding pre-main sequence stars (e.g.. lNatta et al. Il2000l : IScholz et al. II2006I ). we 



198ll). Consistent with observatio ns of 



scale the reference surface density with the stellar mass, So(M^) = 0.18 (M^/Mq) g cm 
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Table 1 lists the ranges in M^ and Xm we consider. The table also lists the main sequence 
lifetime, t^s, defined as the time to reach co re hydrogen exhaustion in the X = 0.71, Y = 
0.27, and Z = 0.02 stellar evolution models of iDemarque et al. I (120041 ). where X is the initial 
mass fraction of hydrogen, Y is the mass fraction of helium, and Z is the metallicity. For 
most of our calculations, the number of annuli in the disk is A^ = 64. To check these results, 
we also calculated models for disks with A^ = 32 around 1 Mq stars. Because the growth of 
planets has large stochastic variations, we repeated the calculations 5-12 times for each set 
of starting conditions, M^, N, Xm, and Qb- Table 1 lists the number of calculations for each 
(M^, Xm) pair. 

Our calculations follow the time evolution of the mass and velocity distributions of 



objects with a range of radii, Tj^ 



to r 



ik 



The upper limit Vmax is always larger 



than the largest object in each annulus. To save computer time in our main calculation, we 
do not consider small objects which do not affect significantly the dynamics and growth of 
larger objects, Tmin = 100 cm. Erosive collisions produce objects with Tj^ < Vmin which are 
'lost' to the model grid. Lost objects are more like ly to be ground down into s maller objects 
than to collide with larger objects in the grid (see lKenyon &: Bromley Il2002al . 2004a). 



To estimate the amount of dusty debris produced by planet formation, we perform a 
second calculation. Each main calculation yields Mj(t), the amount of mass lost to the grid 
per annulus per timestep, and i?io(i), the scale height of the smallest particle (r = rmin) in 
each annulus of the coagulation grid. Objects with sizes smaller than rmin contain a small 
fraction of the mass in each annulus ; thus, the scale height for objects with r < rmin is -f^io(^) 
(JGoldreich. Lithwick. fc Sari II2004I ). The total amount of mass lost from the planetesimal 
grid in each timestep is M{t) = ^^^i Mi(t). The debris has a known size distribu tion, 



nl- 



ntoTi 



where /3 is a constant (see IStern fc Colwell 1 119971 : iKenyon fc Luu I Il999l . and 



references therein). The normalization constant n'-Q depends only on j3 and M{t), which we 
derive at each timestep in the main calculation. To evolve the dust distribution in time, 
we use a simple collision algorithm that includes Poynting-Robertson drag and radiation 
pressure). The optical depth r of the dust follows from integrals over the size distribution in 
each annulus. The optical depth and a radiative transfer solut ion then yield the luminosity 
and radial surface brightness of the dust as a function of time. iKenyon &: Bromley I (l2004al ) 
describe this calculation in more detail. 

Throughout the text, we use simple scaling relations to show how our results depend on 
initial conditions and the properties of the grid. For each set of calculations (Table 1), we 
derive median results for the size distribution, the size of the largest object as a function of a 



^Because the coUisional cascade begins after the gas disk dissipates, we ignore gas drag. 
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and t, and other physical variables. Substituting the inter-quartile range for the dispersion, 
we then perform least-squares fits to relate these median results to input parameters (e.g., 
Xm) and the properties of the grid (e.g., a). For parameters where analytic theory predicts a 
relation (e.g., the growth time as a function of a), we derive the best-fitting coefficient, test 
whether different fitting functions provide better fits to our results, and keep the result that 
minimizes x^ P^r degree of freedom. When analytic theory provides no guidance, we derive 
fitting functions that yield the sensitivity of our results to all important physical variables. 
Thus, our fits test some aspects of analytic theory and guide other aspects. 



3. PLANET FORMATION CALCULATIONS 

3.1. Icy Planet Formation in Disks around 1 M© Stars 

We begin with a discussion of planet formation in disks at 30-150 AU around a 1 Mq 
star. For most disks around low mass stars, the timescale for planet formation is shorter 
than the main sequence lifetime. Thus, the growth of planetesimals into planets and the 
outcome of the collisional cascade depend more on the physics of solid objects than on 
stellar physics. Here, we review the stages in planet growth and describe the outcome of 
the collisional cascade. For a suite of calculations of planet formation in disks of different 
masses, we derive basic relations for the growth time and the maximum planet mass as a 
function of initial disk mass. We also show how the dust production rate and the mass in 
small objects depend on initial disk mass and time. 

The next section compares these results with calculations for 1.5-3 Mq stars. For disks 
around more massive stars, the planet formation timescale is comparable to the main se- 
quence lifetime. Thus, the central star evolves off the main sequence before planet formation 
and the collisional cascade reach a n atural end-state. Durin g post-main sequence evolution, 
the star bri ghtens considerably (e.g ., lDemarque et al. 112004 ) and develops a powerful stellar 



wind (e.g.. iKnapp &: Morris Ill985l . and references therein), melting icy objects in the inner 



disk and ejecting small grains throughout the disk. Thus, we halt our calculations when the 
star evolves off the main sequence. Here, we show how the physics of main sequence stars 
changes the results derived for planet formation around 1 Mq stars. 



3.1.1. Growth of Large Objects 



The formation o f icy planets in the outer regions of a quiescent p lanetesimal disk has 
three distinct stages (JKenyon &: Luu Ill999l : iKenyon fc Bromley Il2004al ). Planetesimals with 
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r < 1 km grow slowly. As they grow, dynamical friction damps e for the l argest objects; dy- 



nami c al friction and viscous stirr i ng raise e for the smallest objects (e.g., iGreenberg et al. 
1984J : IWetherill fc Stewart I Eoosl : IColdreich. Lithwick. fc Sadlbooj ). After ~ 0.1-1 Myr, 
gravitational focusing enhances the collision cross-sections by factors of 10-100. Slow, or- 
derly growth ends. Runaway growth begins. At the inner edge of the disk, the largest objects 
take ~ 3 Myr to grow to ~ 100 km and ~ 30 Myr to grow to ~ 1000 km. Throughout run- 
away growth, the gas disk dissipates. Thus, velocity damping by gas drag ceases; dynamical 
friction and viscous stirring increase e for the smallest objects. Stirring reduces gravitational 
focusing factors, slowing the g rowth of the largest o 



ative to leftover planetesimals (llda fc Makino 



growth ends; oligarchic growth begins (JKokubo fc Ida 



1993 



jjects relative to one anoth er and rel- 
Wetherill fc Stewart I Il993h. Runawav 



1998 



Rafikov 



2003 



Chambers 2006 



Nagasawa et al. 1120071 ). After 30-100 Myr, the largest objects - oligarchs - grow slowly and 
contain an ever increasing fraction of the remaining mass in the disk. 

During the transition from runaway to oligarchic growth, collisions start to produce 
copious amount of dust. Once oligarchs rea ch sizes ~ 500 km, co llisions between 1-10 km 



objects produce debris instead of mergers (JKenyon et al. I l2008l . and references therein). 



Once fragmentation begins, continued stirring leads to a collisional cascade, where leftover 
planetesimals are ground to dust. For dust grains with sizes > 1 nm, the collision t ime 



is much shorter than the tir ne to remove particles by ga s drag (jAdachi et al. I Il976l ) or 



by Poynting- Robertson drag (JBurns. Lamy. fc Soter Ill979l ). Thus, the cascade proceeds to 
partic le sizes S. 1-10 nm, where rad iation pressure removes material on the djTiamical time 
scale ( JBurns. Lamy. &: Soter Ill979l ). Because runaway growth leaves most of the mass in 
1-10 km objects, the collisional cascade effectively removes a significant fraction of the solid 
material in the disk. 

Fig. [Ushows the time evolution of the eccentricity and the mass distributions for objects 
in the innermost 8 annuli of a disk with initial mass distribution similar to the MMSN. To 
minimize stochastic variations, these plots show median results for 15 calculations. During 
slow growth and the early stages of runaway growth, dynamical friction damps e for the 
largest objects and raises e for the smallest objects (Fig. [U right panel, 10 Myr). The mass 
distribution develops a pronounced shoulder from 10 km to ~ 300-500 km. As the evolution 
proceeds, growth concentrates more mass in the largest objects; stirring excites the orbits 
of the smallest objects. After 100 Myr, the collisional cascade removes mass efficiently from 
the smallest objects but leaves the mass distribution of the largest objects unchanged. By 
~ 5-10 Gyr, almost all of the small objects have been removed. 

In these calculations, the rate of planet formation is very sensitive to semimajor axis 
(Fig. [2]). For collisional processes, the growth time in the disk is t oc P/S, where P is the or- 
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bital period (see the Appendix; also lLissauer 1119871 : iKenyon fc Luu lll998l : lGoldreich. Lithwick. fc Sari 



20041 ) ■ For P oc a^/^ and S oc x^a"^/^ (Eq. ([27])), the growth time is t ex a^a;;;^ Thus, 
ahhough it takes only ~ 10 Myr for the largest planets to reach radii of 300-600 km at 
30 AU, the largest objects at a > 100 AU still have r ~ 3-5 km. By 100 Myr, 100 km 
objects form at 75-80 AU. After 1 Gyr, 100 km objects form beyond 125 AU. By the end of 
the calculation at 10 Gyr, 1000 km objects form throughout the disk. 

The formation rate also depends on the initial disk mass (Fig. [3]). For an expected 
growth time t oc a^x^, planets grow faster in more massive disks. At 100 AU, planets with 
r ~ 2000 km form in a massive disk [xm = 3) within 1 Gyr. In a low mass disk with Xm ~ 
1/3, the largest planet at 100 AU grows to r ~ 300 km in 1 Gyr and r ~ 2000 km in 10 
Gyr. For all simulations of disks around 1 Mq stars, the median timescale for the formation 
of the first 1000 km object is 

tiooo ^ 475 x-i-i^ {^^y Myr . (28) 

This relation fits our results for the median growth time to ~ 5% for a = 30-150 AU and 
for Xjn = 1/3 to 3. For each initial disk mass, the inter-quartile range for the formation time 
is ~ 20%. Thus, there is a modest range of outcomes for identical starting conditions. 

In equation [2Sl there is a small difference between the result expected from simple 
theory (t oc x~^) and the result derived in our calculations (t oc x~^'^^). We show in the 
Appendix how gas drag during runaway and oligarchic growth can modify the simple theory 
and explain the result of our calculations. 

Although the timescale to produce the first 1000 km object is a strong function of initial 
disk mass and semimajor axis, the evolution at late times is less sensitive to the starting 
conditions. To derive a simple relation for the median size Vmax of the largest object as a 
function of initial disk mass and semimajor axis, we adopt a simple function 

Tmaccia) = ro e-^^^/''^^"' (29) 



and use an amoeba algorithm (jPress et al. Ill992l ) to derive the fitting parameters Oq, tq, and 



Or as a function of time. For stellar ages 1 Gyr ^ t^ ^ 10 Gyr, the ensemble of calculations 
yields 

ro ^ 1650 xl' [j^] km (30) 

for the radius of the largest object, 

ao ^ 190 xl' [YGfrJ ^^ ^^^^ 
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for the scale length, and a,. ~ 5-6 for the exponent. These relations match our results to 
±5%. The uncertainties are ±3% in the coefficients and ±5% in the exponents. 

These calculati ons produce relatively l ow mass icy planet s with radii ~ 50% larger than 



the radius of Pluto (JYoung fc Binzel Ill994j : lEUiot et al. 



200l 2007). Although these objects 



form relatively rapidly in the inner disk, they grow very slowly at late times. Between 1- 
10 Gyr, most large object grow by ~ 10-20% in radius (~ 50% in mass). Because the size of 
the largest object depends weakly on the initial disk mass, nearly all disks form Pluto-mass 
planets. These objects stir leftover planetesimals effectively; thus, nearly all disks develop a 
coUisional cascade. 

Despite the general formation of Pluto-mass planets in any disk, the lowest mass disks 
{xm ^ 1/3) form Plutos inefficiently. The scale length, oq > 150 AU, derived from our 
calculations exceeds the outer radius of the disk. Thus, planet formation does not proceed 
to completion at large a for the lowest mass disks. In these systems, the largest icy planets at 
a ^ 125-150 AU are factors of 3-10 smaller than tq. The large exponent, a,. ~ 5-6, derived 
in our ffis implies a rapid transition - ~ 10-20 AU - between disk regions where the largest 
objects are planets with r^ax ~ ^o and where the largest objects have r^ax ^ 300-500 km. 

In our calculations, the collisional cascade limits the size of the largest objects. Once 
a few objects have radii > 1000 km, they stir up leftover planetesimals to the disruption 
velocity. When the collisional cascade begins, the timescale for 1 km planetesimals to collide 
and fragment into smaller objects is shorter than the timescale for oligarchs to accrete leftover 
planetesimals. Thus, the growth of the largest objects stalls at ~ 1000-2000 km (~ 0.01- 
0.02 M^). Occasionally, runaway and oligarchic growth produce a very large object with 
r ~ 5000 km (~ 0.1 M^), but these objects form in only ~ 5-10% of the simulations. These 
objects form at random semimajor axes and tend to form in more massive disks. 

The large radial variation in the formation time produces dramatic differences in the 
mass distribution as a function of semimajor axis (Fig. Hj). In the inner disk, rapid growth 
leads to many objects with r > 1000 km (Table [2]). With many large objects stirring the 
leftover planetesimals in the inner disk, the collisional cascade removes most of the mass in 
small objects (Fig. HJ left panel). In the outer disk, slow growth results in a handful of Pluto- 
mass objects. A few large objects cannot stir leftover planetesimals efficiently. Thus, the 
collisional cascade is weak and leaves a substantial amount of mass in 1-10 km planetesimals 
(Fig. m right panel). 

The growth of objects as a function of semimajor axis and time is not sensitive to 
the size of the model grid (Fig. [5]). For two sets of calculations with 32 annuli (cyan and 
magenta points), the median radius of the largest object in each annulus is nearly identical 
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to results for calculations with 64 annuli (black points). The results for Eqs. fl28] - [Ml) are also 
independent of the grid. In principle, long-range stirring from planets at small a can influence 
runaway growth of objects at large a. Our results suggest that icy planet formation at large 
semimajor axes is not influenced by the formation of small icy planets at small semimajor 
axes. 

To conclude this discussion of the formation of large objects in a planetesimal disk, we 
quote several simple relations for the amount of solid material in small and large objects as 
a function of initial disk mass and semimajor axis for the ensemble of calculations around a 
1 Mq star. At 3-10 Gyr, the median fraction of solids remaining in the disk is 

For a MMSN with Xm = 1, the amount of mass remaining in the disk at 3-10 Gyr ranges 
from 9% of the initial mass at 30 AU to roughly 50% of the initial mass at 150 AU. Thus, 
the inner disk is substantially depleted, while the outer disk contains a significant fraction 
of its initial mass. 

For each Xm-, the median fraction of the initial disk mass in 1000 km and larger objects 
is 

Aooo = 0.035 (^) . (33) 

The median fraction of the mass in 100 km and larger objects is roughly 50% larger, 

;,„.0.06(5i^V (34) 



For the ensemble of calculations, the typical inter-quartile range is ~ 0.1 for fs and ~ 20% 
for /loo and /looo- Thus, the mass distributions in our calculations are top-heavy, with more 
mass in 1000+ km objects than in 100-1000 km objects. 

These relations demonstrate that planet formation at 30-150 AU is very inefficient. For 
all disk radii in this range, only ~ 6% or less of the initial population of 1 km objects is 
incorporated into large objects with radii exceeding 100 km. In the inner disk (30-50 AU), the 
coUisional cascade is very efficient at removing leftover 1-10 km objects. Thus, at 3-10 Gyr, 
the large objects contain most of the mass in the inner disk. In the outer disk (100-150 AU), 
the collisional cascade does not have enough time to remove leftover planetesimals. Thus, 
small objects with radii of 1-10 km contain most of the remaining mass at 100-150 AU. 
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3.1.2. Evolution of Dust 

At all semimajor axes, the coUisional cascade converts a large fraction of the initial mass 
in solids into small dust grains. Because oligarchs and leftover planetesimals are unobservable 
with current techniques, dust emission provides the sole observational diagnostic of the 
growth of icy planets at 30-150 AU around other stars. Here, we describe the evolution of 
these small particles and demonstrate that the collisional cascade is observable. 

Two physical processes set the visibility of dust grains in a debris disk. Once significant 
fragmentation begins, collisions gradually grind the fragments to dust. When dust grains 
are small enough, radiative processes remove them from the disk. For disks at 30-150 AU, 
radiation pressure dominates mass loss for t < 1-3 Gyr and removes 65% to 70% of the total 
mass loss. Poynting-Robertson drag removes material at late times and is responsible for 
30% to 35% of the total mass loss. Because the gas density is negligible once the collisional 
cascade begins, gas drag is unimportant. 

To describe our results, we divide dusty debris into 'large grains' with 1 mm ^ t < 1 
m, 'small grains' with 1 /xm ^ t < 1 mm, and 'very small grains' with r < 1 /im. Collisions 
dominate the evolution of large grains at all times. For t < 1-3 Gyr, collisions dominate 
the evolution of small grains; Poynting-Robertson drag then removes grains with radii of 
1-100 /xm on Gyr timescales. Radiation pressure removes very small grains on the local 
dynamical timescale. Thus, radiation pressure produces a 'wind' of very small grains in the 
disk midplane. 

Fig. [H] shows the time evolution of the dust production rate for very small grains as a 
function of initial disk mass. At the start of each calculation, dynamical friction and collisions 
damp orbital eccentricities. Thus, collisions produce less and less debris; the dust production 
rate declines with time. As oligarchs reach radii of ~ 500 km, they stir leftover planetesimals 
along their orbits. Dust production increases. Because oligarchs continue to grow, they stir 
leftover planetesimals to larger and larger velocities. Collision energies rapidly exceed the 
critical disruption energy; the dust production rate increases dramatically (Eq. HHS]). When 
oligarchs start to reach their maximum radii in the inner disk (Eq. (I29l) ). the dust production 
rate peaks. As oligarchs grow at larger and larger disk radii, the dust production rate slowly 
declines. 

Although the outer disk contains more mass, the global dust production rate declines 
with time for two reasons. Large oligarchs form at late times in the outer disk (Fig. [2]), but 
the smaller disk surface density and the longer orbital periods lead to smaller collision rates. 
Smaller collision rates yield smaller dust production rates. Initially, collisions dominate 
Poynting-Robertson drag; thus, radiation pressure ejects the smallest grains on the local 
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orbital timescale (e.g., iKrivov et al. I I2OOOI : IWyatt I l2005l ) . The dust production rate then 
dedines roughly as t~^. At late times, the collision rates decline. Poynting- Robertson drag 
then removes larger grains from the disk, which reduces the population of very small grains 
from ero sive collisions. The dust production rate then declin es with evolution time as t~^ 
see also Dorninik &: Decin 20031: Kenvon fc Bromley 112004a . 



200a IWvatt 



2005 



Wvatt et al. ll2007al Jbll 



The time evolution of the collision rate in the disk yields a simple relation between the 
maximum dust production rate Mmax and the initial disk mass. For the complete set of 
calculations, 

Mma. ~ 6.5 X 10^° xl g yr-1 . (35) 



For a MMSN with Xm = 1, the maximum rate is ~ 0.1 M® every million years. The collision 
rate scales with the square of the number density of objects; thus, the dust production rate 
grows as the square of the initial disk mass (e.g., M^ax oc x^). 

The timescale for the peak in dust production is shorter than the timescale for the 
production of 1000 km objects in the disk. 






14 x-i Myr . 



(36) 



This time corresponds roughly to the time of peak stirring of leftover planetesimals by 
oligarchs growing in the inner disk, starting the coUisional cascade. Oligarchs form faster in 
more massive disks; thus, the dust production rate peaks earlier in more massive disks. 

Fig. [7] shows the time evolution of the mass in small grains for disks with a range of 
initial masses. Initially, the dust production rates are small (Fig. [H]) and the dust mass 
increases slowly with time. Once large oligarchs form in the inner disk, the dust production 
rate - and thus the mass in small grains - grows rapidly. For all disks, it takes ~ 5-10 Myr 
for the mass in small grains to grow from 10~® M® (which is unobservable with current 
techology) to ~ 1 — 10 x 10~^ M® (which is observable). After this rapid rise, oligarchs form 
at larger and larger disk radii, leading to enhanced dust production farther and farther out 
in the disk. The dust mass then rises slowly with time. Once oligarchs form at the outer 
edge of the disk, the collisional cascade removes more and more solid material throughout 
the disk. The dust mass then begins to decline. 

The maximum mass in small grains scales with the initial disk mass. 



M, 



max, small 



0.013 x™ Mff 



(37) 



Because the derived size distributions are dominate d by collisional processes, the maximum 



mass in large grains is roughly 40 times larger (e.g., iDohnanyi Ill969l : IWilliams &: Wetherill 
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1994 : iKenvon fc Bromlev Il2004cl : iKrivov et al. 1 120061 ) . 

^''J''max,large ^~~^ U.O Xm -'^^S 



(38) 



In both cases, the larger colhsion rate in more massive disks leads to more dust. Although 
these dust masses are significant, they are small compared to the mass in objects with r > 
100 km. The typical mass in large grains is < 10% of the mass in 100 km and larger objects. 
The mass in small grains is < 0.25% of the mass in the largest objects. 

The timescale to reach the maximum dust mass is a function of the particle size. For 
the small grains, 

trr,.n...r..nu ~ 270 xz}^^ Myr . (39) 



For the large grains. 



^max,smaU ~ 270 xj^^ Myr . 
tmax,large ~ 600 X;;^/^ Myr . 



(40) 



Several features of the coUisional cascade set these timescales. Early in the evolution, the 
collision timescale for all particle sizes is smaller than the timescale for Poynting-Robertson 
drag. Thus, the cascade erodes small particles until radiation pressure ejects them. As planet 
formation propagates out through the disk, collisions produce more and more small grains. 
Because the mass in grains is set by a balance between the collision time, which scales as x~^ 

— 1/2 

and the local dynamical time, which scales as Xm , the timescale to reach the maximum 

— 1/2 

grain mass is proportional to Xm ■ As the collision rate declines, Poynting-Robertson drag 
starts to remove mass from the disk. This drag removes smaller particles from the disk more 
effectively than it removes large particles. Thus, the mass in small grains peaks before the 
mass in larger grains. 

These results suggest that the mass in collisional debris is large, roughly a lunar mass 
in 0.001-1 mm grains, throug hout the lifetime of a 1 M(^ star. Although the Solar Sys- 
tem has much less dust (e.g., iLandgraf et al. I l2002l : iNesvorny et al. II2006I . and references 
therein) , these large disk masses are compar able to the mass in dust grains detected in 



many debris disks aroun d other stars (e.g., iBeichman et al. I l2006l : iTrilling et al. 1 12008 



Moro-Martin et al. 1 120071 ). Because our dust production rates are observable, the large 



range in dust masses as a function of initial disk mass and time implies a correspondingly 
large range in the observable properties of debris disks, such as the disk luminosity and the 
IR excess, at fixed stellar age. Because the dust production rate declines with time, we ex- 
pect the disk luminosity and IR excesses to decline with time. We derive detailed predictions 
for this evolution in §4 and compare these results with observations in §5. 
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3.2. Icy Planet Formation in 1.5-3 Mq Stars 

Stellar evolution is an important feature of icy planet for mation at 30-150 AU. Because 



the m ain sequence lifetime {tms oc M^ "•, with n ~ 3-3.5; e.g., Ilben 1119671 : iDemarque et al 



20041 ) is more sensitive to stellar mass than the timescale to produce large icy planets [t oc 
M^ ; see below), massive stars evolve off the main sequence before oligarchic growth and 
the coUisional cascade remove solid material in the outer disk. After a 1-3 M0 star evolves 
off the main sequence, it becomes more luminous (as a red giant) and hotter (as a post- 
AGB star). During this evolution, icy planetesimals and planets melt, decreasing collision 
rates and changing the outcome of the collisional cascadcl. Short main sequence lifetimes 
of massive stars thus lead to clear differences in the amount of solid material in large and 
small objects in the outer disk. 

The stellar mass also affects the outcome of icy planet formation. The timescale for 
planet formation scales with orbital period and the surface density, t oc P/E (see the Ap- 
pendix). For a disk with S = Sq x^ a~^/^ (Eq. fl27|) ) and fixed Sq Xm, the formation time is 
t (X a^ Mi, . Thus, at fixed a in disks with identical surface density distributions, planets 
form faster around more massive stars. However, disk masses in the yo ungest stars scale 



with stellar mass (e.g., Md oc M^; iNatta et al. II2OOOI : IScholz et al. II2OO6I ). Thus, Eq scales 



with stellar mass, Sq oc M^ oc M^. Combining these relations leads to a formation time 
t (X a^ Mj7 • Thus, at fixed a in typical disks, icy planets form ~ 5 times faster around 3 
Mq stars than around 1 M0 stars. 

To illustrate how stellar mass and evolution affect planet formation, we begin with the 
growth of large objects at 40 AU and at 100 AU (Fig. [S]). For disks with identical initial 
surface density distributions, planets at the same a in disks around 3 Mq stars grow ~ 1.7 
times faster than planets around 1 Mq stars. Fig. [8] also shows the clear scaling of growth 
time with semimajor axis, t oc a^ for a disk with E oc a~^/^. The simple scaling of the growth 
time with disk mass and orbital period in our calculations leads to a general relation for the 
median timescale for the formation of the first 1000 km object in disks at 30-150 AU, 



2 Mq 



3/2 



For 1-3 Mq stars, this relation fits our results for the median growth time to ~ 7% for a = 
30-150 AU and for Xm = 1/3 to 3. For all initial disk masses, the inter-quartile range for the 



^We assume that melting is accompanied by a loss of volatiles and an increase in the mass density of 
leftover planetesimals. 
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formation time is ~ 20%. Thus, there is a modest range of outcomes for identical starting 
conditions around 1-3 Mq stars. 

Aside from the extra factor of a;~°'^^, this relation follows the prediction of t oc x^M^ ' a^ 
from standard coagulation theory. As outlined in the Appendix, velocity damping from gas 
drag can speed up planet formation in more massive disks. 

Fig. [9] further shows how the growth time varies with a and M^ for constant Xm- At 100 
Myr, icy planets are close to their maximum sizes at 30-50 AU in the inner disk. At large 
disk radii (a ~ 100-150 AU), planet formation is clearly faster around more massive stars. 
For all stars with M^ = 1-3 Mq, the rmax{0') relations have a similar morphology, consisting 
of a plateau at Vmax ~ 1000-2000 km and a steep decline of Vmax with increasing a. As in 
§3.1.1, we fit our results to a simple relation between r^ax, O'l and time (Eq. fl2I?]) ). For all 
of our calculations, we derive an exponent a^ ~ 5-6 and a characteristic maximum radius 

with the exponent b = 0.22 + 0.033 M^/Mq. The disk scale length is 

/ 3t \°-^ 
ao ^ 190 x^M -^ AU (43) 

V ''■ms / 

These results are valid for late times, t ^ 0.1-1 t^s- 

Eq. (142|) shows that the maximum sizes of icy planets at 30-150 AU are relatively 
insensitive to initial disk mass, stellar mass, or time. For disks with identical Xm around 
1-3 Mq stars, the largest icy planets around 3 M© stars are only ~ 10%-20% larger than 
the largest icy planets around solar-type stars. This range is comparable to the range in 
maximum sizes for planets formed in identical disks around stars of identical mass (§3.1.1). 
Disks with a factor of 10 range in Xm yield planets with a 20%-30% range in radii (a factor 
of ~ 2 in mass). Thus, our calculations predict that the largest icy planets at 30-150 AU 
around 1-3 Mq stars have masses comparable to Pluto and other large Kuiper belt objects 
in the Solar System beyond the orbit of Neptune. 

Although the maximum sizes of planets are fairly independent of initial conditions, the 
number of Pluto-mass objects np is sensitive to disk mass and stellar mass (Table [3]). In 
the inner disk (30-60 AU), np scales roughly with initial disk mass and stellar mass. In 
the outer disk (100-150 AU), the formation timescale for icy planets is comparable to the 
main sequence lifetime. Thus, np scales with initial disk mass and stellar mass only for the 
most massive disks. In lower mass disks, stars evolve off the main sequence before disks can 
produce large numbers of Pluto-mass objects. 
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Cumulative mass distributions provide another useful comparison of icy planet formation 
as a function of disk mass and stellar mass. For disks with identical initial surface density 
distributions at a = 30-37 AU, the shape of the mass distribution is fairly independent of 
stellar mass at 100 Myr (Fig. [TDl left panel). Because large icy planets form first in disks 
around more massive stars, disks of fixed age around 3 M© stars have more mass in larger 
planets and are more coUisionally depleted than disks around 1 M© stars. 

For calculations in scaled MMSN (Fig. [TDl right panel), growth is a stronger function of 
stellar mass. As predicted by the simple scaling relations, disks around 3 Mq stars have ~ 5 
times more mass in large objects than 1 M© stars. More mass in large objects produces more 
stirring, enhancing mass loss in the coUisional cascade. Despite large difference in initial disk 
mass, the mass distributions of disks around 1-3 Mq stars are very similar for r < 1 km at 
similar times. 

Although planet formation proceeds faster with increasing stellar mass, stellar evolution 
halts the collisional cascade and the growth of planets in more massive stars (Fig. [TTj) . 
Planets reach their maximum radii in the inner disks for all 1-3 Mq stars; however, the 
timescale for the collisional cascade to run to completion is long compared to the main 
sequence lifetimes of 2-3 Mq stars. Thus, the collisional cascade removes a larger fraction 
of material from the inner disks around 1 Mq stars than from the inner disks of 2-3 Mq 
stars (Fig. [TTl left panel). In the outer disk, the growth time for 1000 km planets is large 
compared to the main sequence lifetime for 2-3 Mq stars. Thus, planets do not reach their 
maximum radii of ~ 1000-2000 km in the outer disks of 2-3 Mq stars (Fig. [TTl right panel) 

To conclude our comparison of icy planet formation around stars of different masses, we 
consider the long-term evolution of all solids in the disk. For t = 0.3 — 1.0 tms, the median 
fraction of solids remaining in the disk is 

1/4 



■' VlOOAU/ \XmJ 



for the ensemble of calculations for 1.5 Mq stars and 



1/3 



for calculations for 2-3 Mq stars. At 30-50 AU, all stars lose roughly the same fraction of 
mass from the disk. At larger disk radii, more massive stars evolve off the main sequence 
before the collisional cascade removes most of the leftover 1-10 km planetesimals. Thus, icy 
planet formation around lower mass stars converts a larger fraction of the initial solid mass 
into dusty debris. 
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Although planet formation around massive stars converts a smaller fraction of the initial 
mass into dusty debris, icy planet formation is equally efficient at producing massive objects 
around all 1-3 Mq stars. For all disks in our calculations, the median fraction of the initial 
disk mass in 1000 km and larger objects is 

Aooo = 0.035 f^^^l . (46) 



The median fraction of the mass in 100 km and larger objects is ~ 50% larger, 

/,„.0.06f?i^V (47) 



Mass distributions for icy planets are top heavy for all 1-3 M0 stars. As for calculations 
around 1 M0 stars, the typical inter-quartile ranges are ~ 0.1 for fs and ~ 20% for /loo and 
/looo- Thus, identical starting conditions lead to a modest range of outcomes. 



3.2.1. Evolution of Dust 

The evolution of dusty debris in disks around 1.5-3 M© stars generally follows the 
evolution for 1 M© stars. As oligarchs form, stirring leads to a collisional cascade that 
converts 10 km and smaller objects into small dust grains. Because planets form more 
rapidly around more massive stars, disks around massive stars produce dust sooner than 
disks around lower mass stars. In our calculations, the disk mass scales with the stellar 
mass. Thus, disks around massive stars also produce more dust than disks around lower 
mass stars. However, massive stars do not live as long as lower mass stars, preventing the 
collisional cascade from removing all of the small objects from the disk (Fig. [TOl) . Over the 
lifetime of the star, disks around lower mass stars form more dust than disks around more 
massive stars. 

To illustrate these points. Fig. [12] shows the time evolution of the dust production rate 
for a scaled MMSN (x^ = 1) around 1 Mq, 2 M©, and 3 Mq stars. During runaway growth, 
destructive collisions are rare; thus, the dust production rate declines with time. As runaway 
growth ends, there are three distinct phases in dust production: (i) an exponential rise when 
runaway objects start to stir leftover planetesimals in the inner disk, (ii) a long plateau as 
oligarchs form farther and farther out in the disk, and (iii) a long decline in dust production 
as the collisional cascade depletes the disk of 1-10 km objects. 

Because planets grow more rapidly around more massive stars, the exponential growth 
in dust production occurs first around more massive stars. The timescale for the onset of 
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dust production also scales inversely with disk mass; thus, more massive disks produce dust 
faster than low mass disks. 

When runaway growth ends and oligarchic growth begins, the dust production rate 
reaches a clear plateau (Fig. [12]). We define the onset of the plateau phase as the time 
of maximum dust productioiij. For our calculations, there is a simple relation between the 
timescale of maximum dust production and the masses of the disk and central star, 

(2 M \^'^ 

At this time, our simulations yield a simple relation between the maximum dust production 
rate and the masses of the disk and the central star, 

M^a. ^ 3.5 X 10^1 xl [yJ^J g yr"' • (49) 

Each of these scaling laws has a simple physical origin. The maximum dust production 
rate, M^ax, depends on the collision rate, the square of the disk mass divided by the orbital 
period. Thus, Mmax oc S^/P oc S^M^ . For disks where the surface density scales with 
stellar mass (Eq. (1271)). S oc x^M*. Thus, Mmax oc x^MJ . The timescale to reach this 
rate is the orbital period divided by the disk mass. Thus, t.v oc S^^M^ ' oc x~}M^ 

Once oligarchs form in the outer disk, the dust production rate declines rapidly. Stellar 
evolution sets the duration of this decline. In massive stars, the short main sequence lifetime 
halts the evolution before the coUisional cascade depletes the disk of 1-10 km objects. Thus, 
the dust production rate declines by roughly an order of magnitude before the central star 
evolves off the main sequence. For lower mass stars, the long main sequence lifetime allows 
the coUisional cascade to remove some material in the outer disk. Thus, the dust production 
rate declines by ~ two orders of magnitude before the central star evolves into a red giant. 
We show in §4 how the long-term evolution of the dust production rate as a function of stellar 
mass produces observable differences in the IR excesses of massive and low mass stars. 

Despite the large differences in dust production rates, there are smaller variations in 
the amount of dust as a function of disk mass and stellar mass. Fig. [13] shows the time 
evolution of the median mass in small grains for scaled MMSN around 1-3 M© stars. Once 
the exponential rise in dust production begins, the dust masses rapidly evolve to similar 



•^Because there are small fluctuations in the dust production rate during the plateau phase, we define the 
maximum dust production as the time when the derivative dM[t)/dt first changes sign. 
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configurations with ~ 10^^ g in small grains. For our set of calculations, the maximum mass 
in small grains is 



M 



max, small 



0.026 Xr 



2 Mq 



Ma 



(50) 



The coefficient in this equation, O.O26M0, is roughly twice the mass of the Moon. The 
maximum mass in large particles is 



M 



max,large 



1.0 x. 



M. 



2 Mp, 



Mffi . 



(51) 



The timescale to reach the maximum dust mass is a function of the particle size. For 
the small grains, 



'^max, small 



For the large grains. 



t 



max^large 



loo X,™ 



300 xZ^^^ 



2 M 



© 



-1 



M, 

2 Mq 



Myr . 



Myr . 



(52) 



(53) 



As described in §3.1.2, the collision rate, the dynamical timescale, and Poynting-Robertson 
drag combine to produce the shorter timescale for smaller dust grains. 



3.3. Limitations of the Calculations 



I n previous papers, we hav e described limitations to multiannulus (IBromlev &: Kenvon 



2006 



1998 



Kenvon fc Bromlev Il200ll . 2002a, 2004b, 2005, 2006) and single annulus (jKenvon fc Luu 



, 1999) coagulation calculations. Here, we review how several of these limitations affect 



results for the simulations described above. 

As long as the statistical assumptions underlying the formalism are met , coagula- 
tion calculations provide a r e asonable represent ation of real collision evolution ( Wetherill 



1980l: iGreenberg et al. 



1984 



Lissauer fc Stewart II 19931: 



1997 



Kenyon fc Luu 



Davis et al. 



19851: barge fc Pellat 



Wetherill fc Stewart 



1998 



1991 



Spaute et al. 



1993l : lstern fc Colwell lll997l : IWeidenschimng et al. 



11991 



Inaba et al. II2OOII ). For calculations at 30-150 AU around 1- 



3 Mq stars, the spacing of mass bins in an annulus and the spacing of annuli in the disk limit 
the accuracy of the results. Our standard mass spacing, S = 2, length ens the evolution time 



by 10% to 20% relative to more accurate calculations with 6 < 1.4 (see lKenyon fc Luu Ill998 
and references therein). The radial resolutio n, Aaj/a,; = 0.025, al s o lengt hens the evolution 



time. Compared to calculations described in iKenyon fc Bromley I (l2004bl ). improvements in 
our treatment of interactions among particles in neighboring annuli reduce lags by a factor 
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of two, from ~ 20% to ~ 10%. Combining the lags for mass spacing and radial resolution, 
our evolution timescales are ~ 20% to 30% longer than the actual evolution times. This lag 
is comparable to the dispersion in timescales derived from multiple calculations with iden- 
tical starting conditions. Thus, improvements in resolution are unlikely to alter our results 
significantly. 



3.3.1. Dynamical interactions 

The coagulation algorithm begins to break down when (i) a few large objects contain 
most of mass in the grid and (ii) the gravity of these objects dominates the stirring. For 
r < 500-1000 km, the largest objects contain a small fraction of the mass in an annulus; 
individual dynamical interactions are much smaller than the Fokker-Planck stirring rates. 
Thus, kinetic theory yields good estimates for collisions and stirring among small objects. 
As objects grow beyond ~ 1000 km, however, both assumptions of our statistical approach 
begin to fail: (i) the collisional cascade removes leftover planetesimals, increasing the fraction 
of mass in the largest objects and (ii) individual interactions among the largest objects in 
an annulus can deflect large objects into neighboring annuli, raising collision and stirring 
rates signiflcantly. With ~ 100-1000 Pluto-mass objects at 30-150 AU (see Tables EHS]) , 
interactions among the largest objects could play a signiflcant role in the late-time evolution 
of our models. 

Dynamical interactions among an ensemble of Pluto-mass planets occur when the radial 
spacing is Aa ~ 2\^rH, where rn is the Hill radius in Eq. (|22ll . For planets with mass 
Mp ~ 6 X 10^^ g and M^ = 1 Mq, r^ = 0.001 a. Thus, dynamical interactions among the 
larg est objects are inevitable when up ^ 50-100 in a region wit h a radial extent Aa/a ~ 
0.2 (JGoldreich. Lithwick. Sz Sari II2004I : iKenyon &: Bromley II2006I ). Many of our calculations 
yield such large numbers of Pluto-mass objects. 

To save computer time, we did not ca lculate the typical long-ter m evolution of hundreds 
of Plutos using our hybrid evolution code (JBromley fc Kenyon II2006I ). However, we can infer 
the o utcome from scal i ng the results of calculations for the formation of the Earth at 1 AU 



Chambers II2OO1I : ISromlev fc Kenvon llJooel : iKenvon fc Bromlev l[2006l l At 0.4-2 AU 



(e.g. 

dynamical evolution of 100-200 lunar mass objects produces several Earth-mass planets 
in 10-30 Myr. The spacing of lunar mass objects in these calculations is ~ 70% of the 
critical spacing Aa ~ 2v3r//, similar to the spacing of Pluto-mass objects at late times in 
our calculations at 30-150 AU. Scaling the evolution times by the ratio of orbital periods 
suggests that 100-200 Pluto mass objects collide to form planets with masses ~ 0.1 M^ on 
1-3 Gyr timescales (see also lLevison fc Stewart ll200ll : lGoldreich. Lithwick. fc Sari II2004I ). 
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This analysis suggests that dynamical interactions between large numbers of Plutos at 
30-150 AU are interesting only for low mass stars. For 2-3 Mq stars, the main sequence 
lifetimes are shorter than the time needed for Plutos to interact, to collide, and to grow 
into Mars-mass planets. For lower mass stars, several test calculations with our hybrid 
code confirm that ensembles of 100-200 Plutos can grow into several Mars-mass planets on 
timescales of 2-5 Gycl. Although this evolution leads to some extra stirring of leftover low 
mass planetesimals, there are only small changes in the dust production rate and the total 
mass in small grains. Thus, dynamical interactions among Plutos have little impact on our 
general results. 



3.3.2. Fragmentation parameters 

Fragmentation is another uncertainty in our calculations. We treat destructive collisions 
with an energy- scaling algorithm that uses (i) the ratio of the center of mass collision energy 
Qc to the critical disruption energy Q*j. and (ii) a simple power-law size distribution to 
apportion ejected material into lower mass bins. Although the energy-scaling algorithm 
yields a reasonable treatment of collisions in coagulation calculations, the disruption energy 
Q*^ sets the strength of the coUisional cascade. Large Q*^ leads to a weak cascade with little 
debris; small Q*^ allows a strong cascade with significant debris. Because Q*^ and the size 
distribution of the ejecta set the amount of material in small grains, we now discuss how our 
choices for these input parameters affect our results. 

Detailed comparisons of various approaches suggest that the size distribution of the 
ejected mass has little impact on our results. For the large collision rates in our calculations, 
all methods for dividing ejected mass among lower mass bins - including dividing the ejected 
mass equally among 2-3 lower mass bins - leads to a power-law mass distribution with a 
characteristic slope of dn/dm ~ —0.8 (JDohnanyi Ill969l : IWilliams &: Wetherill 1119941 ). Thus, 
the adopted mass distribution for the ejecta is relatively unimportant. 



Kenyon fc Bromley I (120051 ) and iKenyon et al. I (120081 ) describe how the form of Q*^ in 



Eq. ([5]) impacts collisional evolution of icy objects at 30-150 AU. Here, we divide Q*^ into 
a bulk component (the first term of Eq. ([5])) and a gravity component (the second term 
of Eq. (IS])). Gravity provides nearly all of the binding energy for large objects with r > 
10 km; the bulk component of Q*^ provides most of the binding energy for small objects with 



''For a > 75 AU, the escape velocity of Mars-mass planets exceeds the orbital velocity. Al- 
though dynamical interactions am ong; Mars-mass objects can produce ejections in these circumstances 



([Goldreich. Lithwick. fc Sari 112004 ). damping by leftover planetesimals limits ejections in our simulations. 
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r < 1-10 km. For icy objects with maximum sizes r^ax ~ 2000 km, stirring never leads to 
orbital motions large enough to disrupt objects with r > 10-20 km. Thus, our choices for 
the gravity component of Q^ have little impact on our results. Although both components 
of the bulk strength - Qb and /5{, - contribute to Q*^, quoted uncertainties in Pb derived from 
theoretical simulations lead to unimportant variations in Q^ as a function of r. Thus, we 
concentrate on Qb- 

To quantify the impact of Qb on our results, we first consider the evolution of the dust 
production rate and the amount of material in large and small grains. During runaway 
growth, the dust production rates for models with Qb < 10^ erg g~^ are 10%-20% larger 
than dust production rates for models with Qb > 10^ erg g^^. At the same time, the total 
mass in large and small grains is ~ 10 times larger for models with small Qb than for models 
with large Qb- In both cases, models with the smallest initial disk mass have the largest 
differences as a function of Qb- During oligarchic growth, these differences disappear. For 
models with Qb = I — 10^ erg g"\ the dispersion in dust production rates near the time 
of maximum dust production is ~ 5% or less for all disks around 1-3 Mq stars. Although 
the dispersion in the total mass in large and small grains is a factor of ~ 3 during the early 
stages of oligarchic growth, the dispersion in dust m asses declines to 10% or less at late times 



when the dust masses reach their maximum values (JKenyon fc Bromley Il2004bl ) 



The time variation in dust production rate and total dust mass as a function of Qb has a 



simple physical origin (see also lKenyon et al. 1120081 ). During runaway and oligarchic growth 



the colhsion energies of small objects sc ale with the mass of the largest objects in the grid (see 



also iGoldreich. Lithwick. &: Sari II2004I ). Thus, small objects have larger and larger collision 
energies Qc at later and later evolution times. Because this feature of the evolution depends 
only on gravitational stirring, it is independent of Qb- Throughout the evolution, the mass 
ejected in a colhsion scales with Qc/Qb (Eq- &)■ Thus, calculations with small Qb eject 
more material at early times compared to models with large Qb, leading to a large dispersion 
in the dust production rate and total dust masses early in the evolution. At late times, all 
calculations produce objects with r^ax ~ 1500-2000 km. These large objects stir all leftover 
small planetesimals to large random velocities, where the collision energies Qc ^ Q^ for all 
Qb- All collisions then lead to copious mass loss, which elir ninates the dispersion in du st 



production rates and total dust masses at late times (see also lKenyon &: Bromley Il2004bl ). 



In addition to the small late-time dispersion in dust production rates and total dust 
masses, our results yield negligible differences in the late-time fractions of mass in large 
objects (r > 100 km) as a function of Qb- The median radius of the largest object and 
the median number of Pluto mass objects are also independent of Qb- Thus, our analysis 
suggests that the fragmentation parameters have a small impact on observable quantities. 
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For low mass disks at t < 10 Myr, destructive collisions between planetesimals with small Qfe 
produce more dust than objects with large Q^- Although these differences are probably large 
enough to be observable, they disappear at late times when planets reach their maximum 
sizes. 



3.3.3. Treatment of small particles 

Our algorithm for deriving the evolution of small particles with rmin ^ 1 m is a final 
uncertainty in our calculations. To follow the evolution of sizes and orbits for large objects in 
a reasonable amount of computer time, we do not calculate the evolution of small particles 
directly. Instead, we use the known production rate of small particles from the detailed 
calculation M(t), an adopted power-law size distribution, and a simple collision algorithm 
to evolve the small particle size distribution with time. 

Because we include radiation pressure and Poynting- Robertson drag in this simple treat- 
ment of coUisional evolution, our predicted size distributions consist of three distinct pieces. 
For particle sizes where the coUisional timescale is shorter than the timescale for Poynting- 
Robertson drag, n oc r~°'^. For very small sizes where radiation pressure ejects grains, we 
adopt n oc r~'^'^ for grains in a constant velocity outflowing wind (see the Appendix). For in- 
termediate sizes, Poynting- Robertson drag can remove grains faster than collisions replenish 
them. Thus, the particle number n — > 0. To conserve mass, we solve a continuity equation 
to derive the number density of grains dominated by Poynting- Robertson drag. 

Although our solution for the evolution of small particles is efficient, it does not con- 
sider how fluctuations in the coll i sion a nd fragmentation rates might modify the size dis- 
tribution. ICampo Bagatin et al. I (Il994j ) note that size- dependent fluctuations can produce 
wavy size distributions for 0.1-10 na m particles. In their simulations of the (3 Pic disk. 



Thebault. Augereau. fc Beust I (120031 ) derive steady-state size distributions with substan- 
tial deficits of 0.1-10 mm particles compared to a standard n oc r~° '^ power-law (see also 



Krivov et al. I l2006l : iThebault fc Augereau 1 120071 : iLohne et al. I l2008l ) . If these deficits are 



typical, then our algorithm seriously overestimates the mass in small dust grains and thus 
the infrared fiuxes of debris disks. 

To check for this possibility, we computed several models with a simple version of our 
multiannulus coagulation code. In these tests, we extracted a complete disk model near the 
peak of the coUisional cascade, extended the lower end of the size distribution from rmin = 
1 m to 1 /im using a power law n oc r~", and continued the calculation for ~ 100 Myr with 
collisions and Poynting- Robert son drag but without our Fokker-Planck velocity evolution. 
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To estimate the range of errors in our simple algorithm, we varied the power law exponent 
for the size distribution, a ~ 0.6-1, the power law exponent for the fragmentation law, 
/9b ~ -0.5-0, and the magnitude of the bulk strength Qi, = 1-10^ erg g~^. For a range of 
disk masses around a 2 M© star, this approach provides a straightforward estimate for the 
accuracy of our results for small particles. 

These tests confirm that the simple collision algorithm yields results reasonably close 
to more detailed coagulation calculations. For models with Pb ~ 0, Qh ^ 10^ erg g~^, 
and a ~ 0.6-1, the derived size distributions are within ~ 20% of those predicted by the 
simple model for all particles with r ^ 0.01-100 mm. Although calculations with Qb ^ 10^ 
erg g~^ yie ld larger deviations fro n a the simple model, these are small compared to those 



quoted by iThebault &: Augereau I (120071 ). Because particles with small (3b are harder to 
fragment, calculations with Pb ^ —0.25 tend to produce smaller departures for a power law 
size distribution than those with Pb ^ —0.25. 

Several features of our calculations combine to minimize wavy size distributions for small 
particles in disks at 30-150 AU. Because icy planet formation is inefficient, the coUisional 
cascade begins when most of the initial disk mass is in 1-10 km planetesimals. Fragmentation 
of the leftovers leads to a very large production rate of 1 m and smaller objects. Continued 
fragmentation o f these objects tends to wash out wavy size distributions produc ed by a 



low mass cutoff (ICampo Bagatin et al. 1 1 1994 IThebault. Augereau. &: Beust I l2003l ) . In our 



Fokker-Planck treatment of velocity evolution, leftover planetesimals are also in dynamical 
equilibrium with larger protoplanets that are 'safe' from fragmentation. Thus, the dust 
production rate from the collisional cascade is well-matched to the dynamical state of the 
system and tends to sustain a power-law size distribution for the smallest objects. 



3.4. Highlights of Icy Planet Formation Around 1—3 Mq Stars 

Starting with a disk of 1 km planetesimals, icy planet formation at 30-150 AU follows 
the same path for all 1-3 Mq stars. This evolution has six main features. 

• It takes 5-30 Myr for runaway growth to produce an ensemble of oligarchs with radii 
of 500-1000 km. Throughout runaway growth, oligarchs stir up the orbits of leftover 
planetesimals. Collisions between leftover planetesimals produce more and more debris. 

• From ~ 10 Myr to the main sequence turnoff, planets slowly grow to a characteristic 
radius. For a broad range of input parameters, the maximum size of an icy planet is 
~ 1750 km at 30-150 AU. Because the timescale for planet formation at 100-150 AU 
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is similar to the main sequence lifetime of a 1-3 Mq star, the inner disk contains more 
1500-2000 km planets than the outer disk. 

As planets grow slowly, a coUisional cascade grinds leftover planetesimals to dust. Early 
on, radiation pressure ejects the smallest grains in an outflowing wind. Later, Poynting- 
Robertson drag also removes larger grains from the disk. In our calculations, radiation 
pressure removes roughly twice as much mass from the disk as Poynting-Robertson 
drag. The timescale for the collisional cascade to remove leftover planetesimals is close 
to the main sequence lifetime of the central star. Thus, the cascade removes more 
material from the inner disk than from the outer disk. 

Icy planet formation is inefficient. In our calculations, icy planets with radii exceeding 
1000 km contain < 3-4% of the initial mass in solid material. Objects with radii ~ 
100-1000 km contain ~ 2-3% of the initial mass. Because short stellar lifetimes limit 
the growth of planets in the outer disk, the mass in large objects declines linearly with 
increasing distance from the central star. Thus, the inner region of the disk contains 
many more Pluto-mass objects than the outer region. 

The dust produced by the collisional cascade is observable. For disks around 1-3 M© 
stars, the maximum mass in small dust grains with radii of 1 fim to 1 mm is ~ 1-2 lunar 
masses. This mass is comparable to the masses derived for the most luminous debris 
disks around A-type and G-type stars. The time evolution of the dust production rate 
and the mass in small dust grains suggest the dust luminosity declines with time. 

Dusty debris is a signature of the formation of a planetary system. This debris is 
present throughout the lifetime of the central star. 



4. DEBRIS DISK EVOLUTION 

To convert our derived size distributions into observable quantities, we perform a radia- 
tive transfer calculation. For each evolution time t, we derive the luminosit y L^ and effective 



temp erature Tj, of the central star from the Y^ stellar evolution models (JDemarque et al. 



20041 ). We then compute the optical depth T{a) of each annulus in our model grid. The 
optical depth allows us to derive the fraction of the stellar luminosity L^, absorbed by each 
annulus. For each grain size in each annulus, we derive an equilibrium grain temperature 
T{r, a) and an emitted spectrum. Summing the emitted spectra over r and a yields the pre- 
dicted spectral energy distribution (SED) and the total dust luminosit y Ld as a function of 



time. The Appendix describes this calculation in more detail (see also iKenyon &: Bromley 



2004bf ) 
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In our calculation of observable quantities, the most important input parameters are 
the smallest stable grain size r2 (also known as the 'blowout' radius; see Backman & Paresce 
1993) and the slope q of the emissivity law for small grains. Although several estimates for the 



minimum grain size suggest r^ ^ 0.5 — 2 Mj^ /xm for 1-3 M^r^ stars (e.g. iBurns. Lamy. &: Soter 



19791 : lArtymowicz I Il988l : iBackman fc Paresce I Il993l : iKimura et al. 1 12002| ) , the coefficient 



and the scaling relation are sensitive to the composition, internal structure, and radiative 
properties of the grains. Because observations allow few tests of this relation, we adopt r2 = 
1 yum for all stars. If more luminous stars have larger r2, our calculations overestimate the 
optical depth in small grains. Thus, we overestimate the mid-IR and sub mm excesses. For the 



emissivity, submm data suggest q k, 0.6-1 from a handful of debris disks (JNajita fc Williams 



20051 : IWilliams &: Andrews 1 12006| ) . To provide some balance for our likely underestimate of 
r2, we adopt q' = 1. Grains with smaller q emit more efficiently at longer wavelengths; our 
models then underestimate mid-IR and submm excesses. 

To describe the evolution of observable quantities with time, we focus on the dust 
luminosity Ld and the excesses at IR and submm wavelengths. The fractional dust luminosity 
Lrf/L^ provides a clear measure of the relative luminosity of the debris disk. For excesses at 
specific wavelengths, we quote the total emission of the disk and the central star relative to 
the emission from the stellar photosphere, Fx/Fxfl. With this definition, disks that produce 
no excess have F\/F\,q = 1; disks where the excess emission is comparable to the emission 
from the central star have Fx/Fxfl = 2. 

We begin this section with a discussion of excess emission for 1 M© stars. After dis- 
cussing results for 1.5-3 Mq stars, we conclude this section with a brief summary. To facilitate 
comparisons of our results with observations. Tables HHH] list results for the fractional dust 
luminosity and excesses at 24-850 /xm. The paper version lists the first five lines of results 
for Xm = 1/3, 1, and 3. The electronic version includes all results for these x^,- 



4.1. Evolution for 1 M© stars 

Fig. [Hlshows the evolution of the fractional dust luminosity Ld/L-^ for disks surrounding 
a 1 Mq star. Early in the evolution, collisions produce mergers instead of debris. For an 
ensemble of growing planetesimals, the surface area per unit mass (and hence the opacity) 
decreases with time. Thus, L^/L^ declines with time. Less massive disks have smaller dust 
masses and smaller dust luminosities. As oligarchic growth begins, the dust luminosity rises 
rapidly and reaches a peak L^/L^ ^ 2 x 10~^ in 30-100 Myr. More massive disks reach 
larger peak luminosities earlier than less massive disks. At late times, all disks converge to 
the same dust luminosity, L^/ L^, ^ 10~^ at ~ 10 Gyr. 
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Despite their small fractional dust luminosities, these disks produce large excesses at 
70 yum (Fig. [151 left panel). For massive disks with Xm = 2-3, the 70 yum excess rises from 
Fjo/FiQfi ~ 2-3 at 3 Myr to Fyq/Fjo^ ~ 30-50 at 30 Myr. Lower mass disks with x^ = 1/3 
to 1/2 produce smaller peak excesses at later times, F-jq/F-jq^ ~ 10 at ~ 100 Myr. For all 
disk masses, the 70 /zm excess is close to its maximum value for a short period when planet 
formation peaks in the inner disk. The excess then declines with time. The rapid decline 
leads to modest excesses at late times, F-jq/Fjq^q ~ 3-5 at ~ 1 Gyr and F-jq/ F^q^ ~ 2 at ~ 
3-10 Gyr. 

The large excesses at 70 /im are a simple consequence of blackbody radiation from 
small grains at 30-50 AU around a solar-type star. With typical temperatures ~ 40-60 K, 
these grains emit most of their radiation at ~ 50-70 jU m. The peak flux from a blackbody 
grain at temperature T is Fx^max oc T^ (I Allen 1 1 19761 ). Setting the total disk luminosity 



Ld oc T^ yields Fx^max oc LdT. At this wavelength, the stellar flux follows a Rayleigh- Jeans 
law, -Fa oc T^A~^. Combining these relations and including a correction factor for inefficient 
radiation from small grains yields a simple relation for the 70 /im flux from the disk and 
central star, 

'Ld/L, 
10- 

For the luminosities in Fig. [TU this relation accounts for the 70 /xm excesses in Fig. [T5] at 
all times. 



i^70/i^70,0 ^ 1 + 10 ( ^^ ) . (54) 



At longer wavelengths, the disks in our calculations achieve larger peak excesses and 
stay close to the peak excess for longer periods of time (Fig. [151 right panel). Disks reach 
their peak excesses at 850 /im on timescales similar to those at 70 /tm, ~ 30 Myr for disks 
with Xm = 2-3 and ~ 100 Myr for disks with Xm = 1/3. The fractional excesses at 850 /im 
are a factor of ~ 2 larger than the excesses at 70 /im. Because the emitting region evolves 
more slowly, these disks are luminous for ~ 1 Gyr and then decline with time. Despite the 
rapid decline, the excesses are signiflcant at late times, with -F85o/-^850,o ~ 3-10 at ~ 3-10 
Gyr. 

The time variation of IR excess also depends on the outer radius of the disk. For solar 
mass stars, grains at 30-50 AU in the inner disk produce most of the flux at 50-100 /im. 
Thus, disks with outer radii of 70 AU and 150 AU produce similar 70 /tm excesses for t < 
30 Myr (Fig. [THl left panel). Once the smaller disk reaches peak emission, the 70 /im excess 
begins a dramatic decline. The larger disk maintains the peak excess for ~ 30 Myr and then 
declines more slowly with time. For t > 100 Myr, the smaller disk is a factor of 2-3 fainter 
at 70 fiia than the larger disk. 

The evolution of disks with different sizes is more dramatic at 850 /im (Fig. [THl right 



-34- 



panel). For typical grain temperatures ~ 20-60 K, long wavelength emission from the disk 
follows the Rayleigh- Jeans tail of a set of blackbodies. The radiation from each disk annulus 
is then oc a?T. Because the outer disk produces more long wavelength emission than the 
inner disk, the 850 yum excess scales with the outer disk radius. For t > 1 Gyr, we derive 
-p85o/-p850,o oc a^^i, where aout is the outer radius of the disk and n ~ 3-4. Thus, doubling 
the outer disk radius increases the predicted 850 /im excess by a factor of ~ 10 at late times. 



4.2. Evolution for 1.5-3 M© stars 

Several factors change the evolution of the dust luminosity and the IR/submm excesses 
in stars more massive than 1 Mf;^. More mas sive stars are hotter; for the Y^ stellar evolution 



isochrones T^ oc Mj, ( iDemarque et al. 1120041 ). Thus, grains in the inner disks around massive 
stars are warmer, emit more short wavelength radiation, and produce bluer colors than disks 
around less massive stars. More massive stars also evolve faster, tms oc M~^. Because the 
evolutionary timescales for solids in the disk are much less sensitive to stellar mass, t oc M~^, 
massive stars have more dust at the end of their main sequence lifetime than low mass stars 
(e.g.. Fig. [T3l) . Thus, these systems have relatively large IR excesses when their central stars 
evolve off the main sequence. 

To compare the evolution of dust emission in debris disks around 1-3 Mq stars, we 
begin with the evolution of the dust luminosity (Fig. [T71) . Planets grow faster around 
more massive stars; thus, the dust luminosity rises earlier for more massive stars. Once the 
collisional cascade begins, the timescale to reach the peak luminosity depends only on the 
initial disk mass and the stellar mass, 

t,,ma. ^ 25 x-^/' (^l^^Myr. (55) 

This timescale is similar to the timescale required to produce the first Pluto-mass object in 
the inner disk (Eq. fj4T|) ). The peak luminosity depends only on the initial disk mass 

Ld,max/ L-k ~ 2 X 10 Xm ■ (56) 

The luminosity remains close to the peak for ~ 10-30 Myr and then declines with time. 

Stellar evolution has a clear impact on the evolution of the dust luminosity. For 1 Mq 
stars, the dust luminosity declines by a factor of ~ 20 before the star evolves off the main 
sequence (Fig. [H]). For 3 Mq stars, the typical decline in Ld/L^, is only a factor of ~ 4. 
Because debris disks have roughly the same peak luminosities, an ensemble of debris disks 
around middle-aged low mass stars should be systematically less luminous than disks around 
middle-aged high mass stars. 
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Stellar physics also produces dramatic differences in the behavior of the 24 /xm excess 
with stellar mass (Fig. [181 lower left panel). At 30-50 AU, the grain temperatures range 
from ~ 40-60 K for 1 Mq stars to ~ 80-120 K for 3 Mq stars. For these temperatures, 
radiation at 24 yum is on the Wien side of the blackbody peak and thus varies exponentially 
with temperature. Our calculations for 1 Mq stars produce very little 24 /zm radiation from 
material at 30-150 AU. However, the peak 24 /im excesses reach -F24/-^24,o ~ 20 for disks 
around 3 M© stars. For all 1-3 Mq stars, our results yield 

log F2A,majF2A,o ~ 0.74 (M, - IMq) + 0.27 {M./Mq) log Xm . (57) 

This maximum flux occurs at roughly the same time as the peak dust luminosity. 

At longer wavelengths, the excesses are less sensitive to stellar mass. Radiation at 70 /im 
is at the blackbody peak for grains in the inner disk. Thus, the inner disk produces most of 
the 70 /im excess. The peak excess is then independent of the stellar luminosity and depends 
only on the total disk mass (Fig. [181 upper left panel), 

i^70,„^a./F70,0 ^ 55 Xr (^) • (58) 

For grains at 30-150 AU, radiation at longer wavelengths is on the Rayleigh- Jeans tail of 
the blackbody. Thus, observations at 160-850 yum probe material throughout the disk. At 
160 /xm, extra emission from hotter grains in disks around more massive stars is balanced 
by more flux from the hotter central star. Thus, the excess is independent of stellar mass 
and depends only on Xm (Fig. [13 upper right panel), 

-F'l60,maa;/-^160,0 ~ 65 X^ . (59) 

At 850 /im, grains in disks around 1 M0 stars are closer to their blackbody peaks than grains 
in disks around more massive stars. Thus, the 850 /im excesses are larger for 1 M© stars 
(Fig. [IHl lower right panel), 

40 xj^9 M^ = 1 Mq 

"850,max 7-^850,0 ~ ^ } (60) 

25 x°;9 M^ = 1.5 - 3.0 M^ 



m 



'0 



At 70-850 /im, the time of peak excess is similar to the maximum in the dust luminosity. 
Thus, all excesses at 24-850 /im peak at ~ 20-30 Myr for 1-3 M© stars. 

Following the peak in the excess at 20-30 Myr, the relative disk luminosity and the 
excesses at 24-850 /im decrease monotonically with time. For this evolu tion, simple debris 



disk models predict a power law decline, L^/L^ oc t " with n ^ 1-2 (e.g. JDominik &: Decin 
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20031 : IWyatt et al. ll2007al Jbl). To compare our results for t > td,max with these predictions, 



we adopt 

/, = L,/L, oc r^-^ (61) 

and 

fx = Fx/Fyo oc t-'^- (62) 

and derive the power law exponents na = d log fd/d log t and n\ = d log fx/d log t from all 
of our calculations as a function of disk mass, stellar mass, and time. 

Throughout the evolution of all our debris disk models, n^ changes continuously with 
time. For t > td^maxi collisions and radiation pressure dominate the removal of small grains. 
As collision rates slowly decline with time, the exponent increases slowly from n^ ^ to 
Hd ^ 1. When the central star approaches the end of its main sequence lifetime, Poynting- 
Robertson drag starts to dominate collisions. The disk luminosity then decreases rapidly; Ud 
increases from ~ 1 to ~ 2. Because most systems are collisionally-dominated, our calculations 
yield a typical nd ~ 0.6-0.8. 

For A ~ 24-850 /im, the exponents nx follow the evolution of n^. Because collision rates 
are larger in the warmer, inner disk than in the colder outer disk, nx increases slowly with 
A. Thus, the typical ^24 ~ 0.6-0.8 is smaller than the typical rigso ~ 0.8-1.0. 

The exponents Ud and nx are somewhat sensitive to the disk mass and the stellar mass. 
At fixed stellar mass, more massive disks evolve faster. Thus, Ud changes faster for more 
massive disks and is larger at the main sequence turnoff. For fixed disk mass, lower mass stars 
live longer and have more time to reach the Poynting-Robertson drag-dominated regime. Our 
results suggest a 0.1-0.2 range in nd and nx for a factor of 10 range in Xm and a factor of 3 
range in stellar mass. 

In addition to excesses at specific wavelengths, the evolution of color excesses yield 
interesting trends with stellar mass and time. Because the 24 /xm excess is sensitive to 
stellar mass, the [24]- [70] color cleanly distinguishes debris disks around stars of different 
masses (Fig. [H]). For 2-3 Mq stars, [24]- [70] rises rapidly to [24]- [70] ^ 1-2 at ~ 1 Myr and 
then rises slowly throughout the main sequence lifetime of the central star. For 1-1.5 M© 
stars, the color rises later, reaches [24]-[70] ~ 3-4 at 30-100 Myr, and then declines slowly. 

For disks at 30-150 AU, the variation of [24]-[70] with M^ depends solely on the prop- 
erties of the central star. Because more massive stars are hotter, their disks are warmer. 
Warmer disks produce bluer colors. Thus, the peak [24]- [70] scales with M^. 

The mass-dependent color evolution of debris disks at 30-150 AU suggests that color- 
color diagrams can discriminate masses of the central star. In Fig. [201 color-color tracks for 
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scaled MMSN around 1.5 M© stars are clearly distinct from tracks for scaled MMSN around 
2 Mq and 3 Mq stars. In Fig. [21], tracks for a range of disks around 2 Mq stars define a 
triangle- shaped locus distinct from the tracks for 1.5 Mq and 3 Mq stars. 

To establish a triangular debris disk locus for each stellar mass, we define two vectors. 
Adopting a vertex, xo,yo, the upper boundary of the locus is a vector connecting the vertex 
with an upper point, x„,y„. The lower boundary is a second vector connecting the vertex 
with a lower point, xi,yi. Table [H] lists our results for the vertex and the upper /lower points 
as a function of stellar mass. For each stellar mass, colors for debris disks at 30-150 AU 
lie within the area defined by the two vectors. More massive disks produce redder colors. 
Within each locus, the initial disk mass scales with distance from the vertex. 

When dust inside ~ 30 AU produces a small IR excess, this color-color diagram provides 
a useful discriminant of stellar mass. For disks at 30-150 AU around 1-3 Mq stars, the typical 
[5.8]-[8] color is small, with [5.8]-[8] < 0.1 at all times. Predicted colors for terrestrial debris 
disks are much large r. For 3 Mq (1.5 Mq) sta rs, we predict maximum colors [5.8]-[8] ~ 



0.5-1 (0.2-0.5) (e.g.. iKenyon fc Bromley Il2004al . 2005). Thus, mid-IR color-color diagrams 



are useful diagnostics of the outer disk for [5. 8] -[8] < 0.1. 



4.3. Summary 

Planet formation and stellar evolution combine to produce several robust trends in the 
time evolution of the dust luminosity and IR/submm excesses from debris disks around 1-3 

Mq stars. 

For scaled MMSN, the maximum dust luminosity is Ld^max ~ 2 x 10^'^. For an ensemble 
of debris disks, the range in the peak dust luminosity scales with the initial mass of solid 
material in the disk. The dust luminosity reaches this peak at roughly the time when the 
first Pluto-mass objects form at 30-50 AU. Following this peak, the luminosity declines as 
t~"d with Ud ~ 0.6-0.8. Because lower mass stars have longer main sequence lifetimes, debris 
disks around lower mass stars reach smaller fractional dust luminosities at late times. 

The IR/submm excesses from debris disks at 30-150 AU are sensitive to the mass of 
the central star. At 24 /im, disks around more massive stars produce larger excesses; disks 
around stars with M^ < 1 Mq produce negligible excesses at 24 /im. At 70 /im, the excess is 
a simple function of the total dust luminosity, F^q/Fjoo ~ 1 + lO^L^/L^. At 850 /xm, debris 
disks around 1 Mq stars produce larger peak excesses than disks around more massive stars. 
At late times, however, the typical 850 fim excess is fairly independent of stellar mass, with 
-^85o/-^850,o ~ 3-5 for stars with ages t r^ t^s- 
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Among stars with different masses, mid-IR colors provide a sensitive discriminant of 
debris disk evolution when the [5.8]-[8] color is small (Fig. [T9142T]) . For 2-3 Mq stars, [24]- 
[70] slowly becomes redder with the age of the central star; for 1-1.5 Mq stars, [24]-[70] rises 
more rapidly, remains at peak color for 300 Myr to 1 Gyr, and then declines rapidly with 
time. For all stars, [8]-[24] and [24]-[70] correlate with stellar mass. Debris disks around 
2~3 Mq (1-2 Mq) stars have redder (bluer) [8]-[24] colors and bluer (redder) [24]-[70] colors. 
Thus, an [8]-[24] vs [24]-[70] color-color diagram provides a way to analyze debris disks 
around stars with different masses (Table [2]) • 



5. APPLICATIONS 



To test whether our predictions provide a reasonable match to observations, we now 
consider several applications of our models to real systems. For these calculations, the 
broad trends in the evolution of IR excesses and colors are sensitive to the physics of planet 
formation and the collisional cascade. Thus, our main goal is to compare our results with 
observed trends of excesses and colors for large samples of main sequence stars observed with 
the IRAS, ISO, and Spitzer satellites. In addition to long-term trends, the absolute level of 
the excesses depends on r2 and q. Thus, our second goal is to learn whether our assumptions 
yield mid-IR and submm excesses similar to those observed. 

We begin with an analysis of Spitzer data for the prototypical debris disk, Vega. After 
demonstrating that our models can explain the mid-IR fluxes and morphology of Spitzer 
images for this system, we show that our predicti ons provide a good match to obser vations 
of mid-IR excesses for a s ample of A-type stars (IRieke et al. II2005I: ISu et al. II2006I ) and a 
sample of solar-type stars ( Beichman et al. II2006I : iHillenbrand et al. II2008I ). 



5.1. The Vega Disk 

Observations of Vega with IRAS first revealed a large excess of emis sion above the A- 
type photosphere for wavelengths exceeding 12 /xm (jAumann et al. Ill984l ). The best-fitting 
single temperature blackbody to the IRAS data yields a temperature of ~ 85 K, a fractional 
luminosity of ~ 2.5 x 10~^ relative to the central star, and a radius of ~ 150-200 AU for 
the emitting material. Because the lifet i me fo r small grains at 150-200 AU is much shorter 
than the age of Vega, lAumann et al. I (Il984l ) concluded that the grains have sizes larger 



than 1 mm. Thus, Vega provided the first direct evidence for grain growth outside the Solar 
System. 
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Since the lAumann et al. I (Il984[) discovery, Vega has become the protot ypical debris 



disk (e.g.. lBacknian &: Paresce Ill993l : lArtymowicz 1119971 : iLagrange et al. Il2000l ). The debris 



consists of a bright torus with small-sca l e clumps at 80- 1 000 AU from t he central star 



f lHoUand et al. I Il998l : IWilner et al. I l2002l : iLiu et al 



2004 



Su et al. II2OO5I ) and a smaller 



20061). Dust in the small disk 



disk of debris at ~ 1 AU from the central star (lAbsil et al. 

is hot (~ 1500 K), luminous [LaLLjl. 

with a diameter of ~ 0.5-1 AU (jAbsil et al. 1120061 ). This dust might be a result of collisions 



5_x 10_^), and mostly confined to a narrow ring 



between larger objects at 1 AU or grains lost from icy comets at 80-100 AU in the outer 
disk. 



Recently, ISu et al. I (120051 ) analyzed high quality Spitzer images at 24, 70, and 160 /im. 
Their results demonstrate that the large-scale debris consists of a bright ring at 80-200 AU 
and a smooth 'halo' that extends to ~ 1000 AU at 160 /im. The halo has an a~^ radial 
density profile, consistent with a wind of small grains ejected by radiation pressure. Fits 
to the radial surface brightness profiles and the spectral energy distribution suggest the 
grains in the wind have sizes of 1-50 /im and a total mass of Mrf^i_5o ~ 3 x 10~^ M^. The 
grains in the bright ring are larger, with typical sizes o f ~ 240 /im, and have a total mass 
of Mrf24o ~ 2 X 10~^ Mq (see also iMarsh et al. II2OO6I ). For an adopted residence time of 
~ 10^ yr in the wind, the mass in small grains implies that larger grains in the ring produce 
smaller dust particles at a rate of ~ 10^^ g s~^. 

To check whether our model predictions can match the MIPS data for Vega, we make a 
simple comparison with median res ults from several calculations. In addition to the observed 
fluxes at 24, 70, 160, and 850 /xm ( ISu et al. 1120051). we adopt publishe d values for the age 
f200 Mvr: ISu et al Il2006h . luminosity (37 L^ : lAufdenberg et al. l [2006[). and mass (2.3 M©; 
Aufdenberg et al. I l2006l ) of the central star. ISu et al. I (120051 ) separate the observed fluxes 
into contributions from the debris disk and the central star. The first row of Table [TU] lists 
these results, along with their derived values for the mass in 1-50 /im dust grains and the 
dust production rate. The rest of Table [TO] lists predictions for four of our debris disk models 
around 2-2.5 M© stars. 

The comparison in Table [10] suggests a reasonable match to the data. Predictions for 
the total mass in 1-50 /im dust grains and the fluxes at 160 /im and at 850 /im bracket the 
observed values. Our models also predict a bright ring in the 70-160 /im dust emission at 
80-130 AU, close to the ob served position of the bright ring inferred from Spitzer images 
(85-200 AU: ISuet aLll2005h . The mass of 0.1-1 mm particles in this ring, ~ 3-5 x lO'^ M®, 
also agrees with the mass in 240 /im grains derived from the Spitzer data. However, our 
models overpredict the fluxes at 24 /im and at 70 /im by a factor of 2-10 and underpredict 
the dust production rate by a similar factor. 
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To understand possible origins for the mismatches between the data and the models, 
we consider the evolution of small grains in our calculations. When the coUisional cascade 
starts removing material from the disk, the most destructive collisions involve grains with 
comparable masses. These collisions gradually erode the parent objects and produce modest 
amounts of debris in smaller particles. Because (i) the collision timescale is much shorter 
than the timescale for Poynting-Robertson drag and (ii) the ratio of the radiation force to 



rpA 



_£xadUU 



grni! 



oc r 



erosion continues until particles reach a 



the gravitational force is /3, 

size r2 where I3rad ^ 0.5-1 (JBurns. Lamy. fc Soter Ill979l ). Particles with r < r2 are ejected. 

For simplicity, we assume r2 ~ 1 /im for all of our calculations. 



In this picture for the coUisional cascade, the ISu et al. I (120051 ) results provide a simple 
solution for the overprediction of the 24-70 /xm fluxes in our calculations. If grains with 
r2 ^ 1 /im are in the wind, our 'Vega models' underestimate the mass in the wind and 
overestimate the mass left behind in the disk. More mass at larger distances from the 
central star lowers the optical depth in the inner disk, reducing the predicted fluxes at short 
wavelengths. Thus, increasing our adopted r2 for Vega models should provide a better match 
between observed and predi cted fluxes at 24 /xm and at 70 /xm. For an adopted L^ = 60 L©, 
Backman fc Paresce I (Il993l ) estimated r2 = 14 /xm. Scaling this result for our adopted L^ = 
37 Lq, r2 ~ 8.5 /xm. Several test calculations with r2 = 10 /im yield predicted 24 /im and 
70 /im fluxes close to the observed values. 

Reconciling the estimated dust production rate with our predictions requires a more 
rigorous analysis of dust producti on and ejection in the Vega debris disk. To derive the dust 
production rate, ISu et al. I (120051 ) assume that (i) the 240 /im grains in the bright ring are 
bound, with Prad ~ and residence times ^ 10^ yr, and (ii) the 1-50 /im grains in the wind 
are unbound, with P^-ad > 1 and residence times ~ 10^ yr. If the coUisional cascade proceeds 
as a gradual erosion of larger objects into smaller objects, however, w e expect a more gradual 
trans i tion from grains wi th jSrad ~ to grains with P^-ad ~ 1 (see also iBurns. Lamy. fc Soter 



19791 : lArtymowicz Ill988l ). Allowing the residence time to change gradually from the bound 
240 /im grains to the unbound 1 /im grains provides a way to lower the apparent dust 
production rate and to resolve the mismatch between our models and the observations. 

To provide an alternate estimate for the residence time of grains in the Vega disk, we 
consider the collision times in the ring and the wind. We adopt the dust masses derived 
from the Spitzer images {Md,i-5o and Md,24o) and a typical particle size (r). For a ring at 
a ~ 150 AU with a width Aa ~ 50 AU, the collision time for a single grain is 



tc~ 



10^ P 



(r) 



10 /im 



(63) 



where P is the local orbital period in yr (JLissauer 1119871 : 1 Wetherill fc Stewart lll993l : lKenyon fc Luu 
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19981 ). Collision times in the wind are similar. For P ~ 10'^ yr at 150 AU, the collision times 
range from ~ 10^ yr for 1 /im grains to 2 — 3 x 10® yr for 240 fim grains. 

If we assume that the residence times are comparable to the collision times, we can 
construct a self-consistent picture for the collisional cascade in the Vega disk. Collisions in 
the bright ring gradually erode 200-300 ^m grains until they reach sizes ^100 /im, when 
they become incorporated into the wind. Collisions in the wind gradually erode the smaller 
grains until they reach sizes ~ 1 /im, when they are ejected rapidly from the system. 

As long as the 200-300 /im grains are replenished from a reservoir of larger grains, this 
cascade can remain in a quasi-steady state over the main sequence lifetime of Vega. The 
required mass for the reservoir of larger objects is ~ 100-1000 times the current mass in 
240 /im grains, ~ 1-5 M^. This mass is small compared to the initial mass of solid material 
in a torus at 80-200 AU in a scaled MMSN, ~ 50-100 M® (Eq. ([27])). Because the optical 
depth of this reservoir is small, it produces a small IR excess compared to the emission from 
smaller grains. 

This picture relies on two features of the collisional cascade. We need an approximate 
equi valence in inass between the large grains in the ring and the small grains in the wind. 
The ISu et al. I (120051 ) mass estimates support this feature. We also need a gradual change 
in grain lifetime from the ~ 10® yr collision timescale of the large grains to the 10^ — 10^ yr 
dynamical lifetime of the smallest grains. Otherwise, collisions in the broad torus cannot 
occur fast enough to maintain the current smooth structure of the wi nd for timescales longer 
than ~ 10^ — 10^ yr. Current theoretical analyses support this idea (IBurns. Lamv. fc Soter 



1979l : lArtymowicz lll988uTakeuchi fc Artymowicz 11200 ItlGrigorieva et al. 1120071 ). Numerical 
simulations of a collisional cascade with a careful treatment of the interactions between the 



radiation field and the small grains could test this proposal in detail (e.g.. lGrigorieva et al. 



2003). 



We conclude that our calculations provide a reasonable match to observations of Vega. 
If we adopt r2 ~ 10 /im, the data are consistent with a standard collisional cascade within a 
broad torus at 80-200 AU. The cascade feeds an outflowing wind of small grains with sizes 
1-50 /im. If the grain lifetime changes smoothly from ~ 10® yr for large grains to ~ 10^ — 10^ 
yr for small grains, the cascade can maintain the wind indefinitely. 



5.2. Debris Disks Around A-type Stars 

Since the discovery of the Vega debris disk, IR AS, ISO, and Spitzer observations have re- 
vealed debris around dozens of nearby A-type stars ( Backman fc Paresce lll993l : lLagrange et al. 
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20001 : iRieke et al. 1 120051 : iSu et al. 1 120061 ). Like Vega, several of these disks are resolved 
and thus provide important information on the radial structure of the dusty disk (e.g., 



Smith fc Terrile I Il984l : IStapelfeldt et al. 1 12004 : iKalas I l2005l : iMever et al. 1 120071 : ISu et al 



20081 ). Although most A-type stars with debris disks are unresolved, the sample is large 
enough to probe the time evolution of debris around 1.5-3 M0 stars. We now consider 
whether our calculations can explain this evolution. 

To compare our model p r edicti ons with observations , we examine data for nearby A- 
type stars from iRieke et al. I (120051 ) and ISu et al. I (120061 ). iRieke et al. I (120051 ) combined 
24-25 /xm data from IRAS and ISO with new 24 fun photometr y from Spi t zer to investigate 



the decay of planetary debris disks around 266 A-type stars. ISu et al. I (120061 ) analyze a 
sample of ~ 160 A-type stars with high quahty 24 /xm and/or 70 /xm data acquired with 
MIPS on Spitzer. The combined sample has 319 (160) sta rs with 24 nm (70 jUm) o bserv ations, 
spectral types B7-A6, and ages 5-850 Myr. From the iKenvon fc Hartmann I (119951) table 



of stellar effective temperatures and spectral types a nd the iDemaraue et al. I (1200411 stellar 
evolution tracks, ~ 75% (~ 85%) of the stars in the iRieke et al. I (120051 ) (ISu et al. II2OO6I ) 



sample have masses of 1.7-2.5 M0. Thus, we compare these data with our results for debris 
disk evolution around 2 Mq stars. 

The observed 24-70 yum excesses of A-type stars show a clear trend with the age of the 
star (Fig. [22] - [23l) . Although the statistics are poor, the data suggest a rise in the 24 /xm 
excess at 5-10 Myr. The larger sample of young stars with 70 /xm excesses provides better 
evidence for this rise. At both wavelengths, the excess has a broad peak for stars with ages 
of 10-30 Myr. At l ater times, the excess declines with time as t~" with n ~ 0.5-1 (see also 



Decin et al. II2003I : iGreaves fc Wvatt II2003I : iRieke et al. II2005I : iRhee et al. Il2007ar ). 



To improve the statistics for 24 /zm excesses around younger star s, ICurrie et al. I (l2008al ) 
added Spitzer data for many young clusters to the lRieke et al. I (120051 ) sample. This expanded 
set of data provides unambiguous evidence for a rise in the typical 24 fim excess at stellar 
ages of 5-10 Myr and a robust peak in the excess at stellar ages of 10-15 Myr. As in Fig. 
22l the 24 fim excesses for this larger sample of A-type stars decline with age from ~ 20 Myr 
to 1 Gyr. 

In addition to the long-term time evolution of mi d-IR excess, the data also indicate a 



large range in the 24-70 /xm excess at fixed stellar age (IRieke et al. II2005I : ICarpenter et al. 



20061 : ISu et al. II2OO6I : ICurrie et al. Il2008al ). Although younger stars are more likely to have 
mid-IR excesses than older stars, there are many stars without excesses at every age. For 
ages < 200 Myr, stars are equally likely to have any excess between zero and the maximum 
excess at that age. As stars age, they are less likely to have an excess close to the maximum 
excess at that age. Thus, the dispersion in the excess declines with time. 



-43- 



Our calculations provide a good match to the time evolution of the amplitude of the 
24-70 fim excesses. At both wavelengths, the models explain the rise in the amplitude at 
5-10 Myr, the maximum at 10-20 Myr, and the slope of the power-law decline at late times. 
For models with x^ = 1^3, the predicted excesses also agree with the maximum observed 
excesses. Although there are a few stars with excesses larger than the model predictions, 
more than 99% of the A stars in this sample have excesses within the range predicted in our 
calculations. 

Our calculations also provide a natural explanation for a large range in the observed 24- 
70 fim excesses at fixed stellar age. At 70 fun, the maximum excess is roughly proportional 
to the initial disk mass (Eq. (I58l) ). Thus, a factor of ten range in initial disk masses yields 
nearly a factor of ten range in the maximum excess at 70 /im. For stars with ages 10- 
300 Myr, the Spitzer observations suggest a factor of ~ 100 range in the 70 /im excess. If 
this range is set by the initial disk mass, our models suggest initial disk masses with Xm = 
0.03-3. 

Variations in the initial disk radius can also produce a range in 24-70 /im excesses at 
fixed stellar age (e.g.. Fig. (TH]). For 2 Mq stars with ages ~ 400-800 Myr, our results suggest 
that a factor of three variation in the outer disk radius (e.g., 50-150 AU) yields a factor of 
two (five) variation in the amplitude of the 24 yum (70 /im) excess. Although the observed 
range in the amplitude of the 24 /im excess for older A stars agrees with this prediction, the 
range at 70 /im is much larger. Thus, variations in the initial disk radius can explain some 
of the observed range of excesses at 24-70 /im. 

Observations of the youngest stars support a large range in initial disk masses and disk 
radii. Submillimeter observations of dusty disks in the nearby Ophiuchus and Taurus- Auriga 
star- forming regions indicate a 2-3 (~ 1) order of magnitude r ange in the masses (radii) of 



disks surrounding young stars with typical ages of ~ 1 Myr (e.g.. lOsterloh fc Beckwith Ill995 



Motte fc Andrei I2OOII : [Andrews fc Wilhams II2005I . 2007a,b). Our models with x^ = 3 have 



disk masses a little smaller than the maximum dust masses derived from the submillimeter 
surveys. Thus, the submillimeter data imply disks with initial masses 0.01 < Xm ^ 5 and 
initial disk radii 50 AU < aout ^ 1000 AU. Disks with this range of initial masses and outer 
radii can produce the range of 24-70 /im excesses observed around nearby A-type stars. 

Variations in the initial surface density distribution can also lead to a range in the 24- 
70 /xm excess. In our calculations, we adopted a 'standard' surface density relation with E oc 
^-3/2^ Compared to this model, disks with shallower (steeper) surface density distributions 
have relatively more (less) mass at large semimajor axes. The outer disk has cooler grains 
than the inner disk; thus, disks with shallower (steeper) surface density distributions should 
produce more (less) flux at longer wavelengths than our standard models. Because the 
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relative fluxes at 24 /iin and at 70 /iin provide a measure of tlie relative disk masses at 
different semimajor axes, color indices provide a natural measure of the gradient of the 
surface density distribution. 

Fig. [2H compares the predicted color evolution for disks around 2 M© stars with data 



from ISu et al. I (120061 ) . Although current observations do not probe the evolution well at 
1-10 Myr, disks have a large color range, [24]-[70] ~ 1-3, for stars with ages ~ 10 Myr. For 
older stars, the data suggest a slow rise in the maximum color from [24]-[70] ~ 3 at 10 Myr 
to [24]-[70] ~ 3.5 at ~ 100 Myr. After ~ 300 Myr, the maximum color declines. For all 
stars older than ~ 100 Myr, the range in color is ~ 3-4 mag. 

Our models match the observed color evolution. At all ages, the predicted colors for 
calculations with Xm = 3 provide a clear upper envelope to the observed colors. The predicted 
colors also explain the slow rise in the maximum observed color for stars with ages of 10- 
100 Myr. To explain the full range in observed colors for 100 Myr to 1 Gyr old stars, we 
require disks with initial masses Xm ~ 0.01-3. This range is similar to the range required 
for the time evolution of the 24 /im and 70 /im excesses. 

The good match to the color observations suggests that the typical initial surface density 
distribution is reasonably close to our adopted S oc a~^/^. For disks with shallower gradi- 
ents, we expect redder colors at later times. A few stars lie above our model predictions; 
however, most stars have bluer colors than models with Xm = 3. Thus, few disks in these 
samples require shallower surface density distributions. Disks with steeper surface density 
distributions can produce stars with blue colors, [24]-[70] ~ 1, at late times. A large sample 
of A-type stars with [70]-[160] colors and spatially resolved observations of the radial dust 
distributions of these stars would provide a constraint on the initial surface density gradient. 

Despite our success in matching these observations, other physical processes may be 
needed to explain the full diversity of debris disk properties for A stars wit h similar ages and 



luminosities. In their analysis of the large debris disk surrounding 7 Oph, ISu et al. I (120081 ) 
examine a dozen main sequence stars with A0-A3 spectral types, ages of 1 50-400 Myr, and 
fractional disk luminosities L^/L^ f» 10~^ to 10~^ (see also lSu et al. II2006I ). Although all of 



these stars have dust with T fa 50-100 K, Fomalhaut has a bright torus of dust with weak 
or negligible emission from a wind of small grains, Vega has a bright torus with a luminous 
wind of small grains, and 7 Oph has an extended disk {aout ~ 500 AU) of dust apparently 
bound to the star. Some A0-A3 stars have warm inner disks with dust temperat ures 



100 K to > 200 K; other A stars have no obvious warm dust emission. ISu et al. I (120081 ) 
conclude that coUisional cascades in disks with a range of masses and other processes, such 
as the formation of giant planets or recent catastrophic collisions, combine to produce the 
wide range of observed properties in this sample. 
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In principle, our models can explain some of this diversity. The ISu et al. I (120081 ) sample 
contains A stars with a factor of five range in L^. Thus, these stars probably have a factor 



offiye range in the blowout radius r2 (see also §5.1; lArtymowicz Ill988l : iBackman &: Paresce 



19931 ). If the proto stellar disks around these st ars had properties similar to those observed 
in Taurus- Auriga (JAndrews fc Williams I l2005l . 2007a), they probably had a factor of ten 
range in initial disk mass, a factor of three range in initial disk radius, and a 50% range 
in the slope of the initial surface density distribution. Co upled with a similar dispe rsion 
in initial conditions for the terrestrial zones of these stars (JKenyon fc Bromley II2005I ). our 
results suggest that this range in initial conditions can produce a broad diversity of debris 
disks. We have started a suite of calculations to address this issue. Larger samples of A 
stars with resolved disks will provide crucial tests of these calculations. 

Other aspects of planet formation are also important. If the cores of gas giant planets 
form before their parent st ars reach the main sequen ce, we expect gas giants at 20-30 AU 
around 2-3 M0 stars (e.g.. [Kennedy fc Kenyon 1120081 ). Gas giants rapi dly remove debris in 



the inner disk and impose struc ture in the debris beyond 30 AU (e.g.. IWilner et al. 112002 



(jJohnson et al. 



Moro-Martm fc Malhotra II2005I ). Because gas giants are common around evolved A stars 



200 



71), gas giants probably play a significant role in the evolution of debris 



disks around A stars. 

Catastrophic col l isions may also produce diversity among A star debris disks (e.g.. 



Wyatt fc Dent I l2002l : ISu et al. I l2005l ). Although debris from complete disruption of co l- 



liding planetesimals is unobservable in our simulations (see also iKenyon fc Bromley II2005I ). 
dynamical events similar to those that produced the L ate Heavy Bombardment in the Solar 
System probably are visible (e.g., IComes et al. II2005I ). Testing this idea requires numeri- 
cal calculations that link the dynamics of rnassiy e planets with the collision al evolution of 
smaller objects (e.g.. ICharnoz fc Morbidelli II2003I : iKenyon fc Bromley 1120061 ). 



We conclude that our debris disk models can explain the overall time evolution of the 
IR excesses and IR colors of A-type main sequence stars at 24 yum and at 70 /im. Our 
calculations for disks with Xm = 1/3 to 3 around 2 Mq stars fit the overall level of the 
excesses and the trends with stellar age. Explaining the full range of observed IR excesses 
and IR colors requires a set of disks with Xm = 0.01-3, as suggested from observations of 
disks around the youngest stars. Matching other properties of these stars - including the 
relative amount of emission from a warm inner disk, an outflowing wind of small grains, and 
a large outer disk - requires calculations that include a broader ra nge of initial disk radii and 



gas giant and terrestr i al planet formation at a,; < 30 AU (e . g.. iKenvon fc Brom 



Nagasawa et al. 1120071 : iKennedv fc Kenvon II2OO8I : llda fc Lin II2OO8I 



Kretke et al. 
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5.3. Debris Disks Around Solar-type Stars 

Although most of the debris disks discovered with IRAS and /5'Q have A- type central 
stars, a few have F-type or G-type central stars with masses of 1-1. 5 Mrn ([Backnian &: Paresce 



19931 : iLagrange et al. II2OOOI : iDecin et al. II2003I : ISong et al. II2005I : iRhee et al. Il2007ar). More 



recent Spitzer observations reveal debris disks aro u nd many solar-type s t ars (IBrvden et al 



2006 : iBeichman et al. I I2OO6I : iMever et al. I I2OO6I : iTrilling et al. I I2OO8I : iHillenbrand et al. 



20081 ). Although several Spitzer programs concentrate on older solar- type stars as prepa- 
ration for detailed planet searches, the range of ages is large enough to provide an initial test 
of our predictions. 

To compare our model predic t ions w ith o bservations, we cons i der da ta for nearby solar- 



type s tars from iBeichman et al. I (120061 ) and IHillenbrand et al. I ( 120081 ) . iBeichman et al. 
(|20061) observed ~ 80 so lar-type stars at 24 /xm and at 70 fim. using MIPS on Spitzer. 
Hillenbrand et al. I (120081 ) analyze ~ 30 stars with 70 /xm excesses out of a sample of 328 stars 



from the Spitzer Leg acy Science Program, "Formation and Evolution of Pl anetary Systems' 



( Meyer et al. II2006I ). After eliminating K-type and M-type stars from the IBeichman et al. 
( I2OO6I ) study, the two programs contain ~ 80 stars with ages of ~ 10 Myr to ~ 10 Gyr. Most 



of these stars have masses of 0.8-1.5 Mq. Thus, we compare these data with our results for 
debris disk evolution around 1 M0 stars. 

The observed 70-160 fun excesses of solar-type stars show trends similar to those ob- 
served in the evolution of A- type stars at 24-70 /xm (Fig. [25] - l26l) . Although the statistics 
for solar- type stars are poor for the youngest stars, the data suggest a rise in the 70-160 fim 
excess at 10-100 Myr. The maximum in the 70 /xm excess is comparable in magnitude but a 
factor of ~ 10 later in time than the maximum 70 /im excess for A-type stars. At 70 fim and 
at 160 /im, the excess follows a roughly power-law decline with time for older stars. Solar- 
type stars also have a large range in excess at all ages, with Fjq/ F-jq^ ^ 1-300 at 100 Myr 
and Fiq/ FjQ^Q ps 1-30 at 1-3 Gyr. For stars with similar ages, the range in the 70 /im excess 
is larger for solar-type stars than for A-type stars. 

Our models match the observed trends for the IR excesses of solar-type stars. At 70 /im, 
65%- 75% of the observations lie within model predictions; at 160 /xm, more than half of the 
observations are within model predictions. For a range of initial disk masses {xm ~ 0.01-3) 
and outer radii {aout ~ 70-150 AU), we predict a large range of excesses at all ages, as 
observed. These model also explain the larger 70 /xm excesses observed for solar-type stars 
relative to A-type stars, the apparent maximum in the 70-160 /im excess at 30-100 Myr, 
and the general power-law decline in the excess flux for the oldest stars. 



Despite this general success, however, the models underpredict the largest observed 
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excesses. At 70 fim, the brightest disks are a factor of 5-10 brighter than disks with Xm 
= 3. At 160 yum, the brightest systems are 3-5 times brighter than our most luminous 
disks. Although the sample of 160 yum sources is small, our models underpredict the largest 
observed fluxes at all ages. 

Changing two assumptions in our mo dels yields a better mat c h to t he observed fluxes at 



70-160 yum. For realistic grain properties, iBurns. Lamy. fc Soter I (119791 ) show that radiation 



pressure from the Sun cannot eject small grains from the Solar System. Reducing the 
minimum stable grain size from r2 = 1 /im to r2 = 0.1 /zm increases our predicted 70 fira 
(160 /xm) fluxes by a factor of 2-3 (1.5-2). Submm observations of several debris disks imply 
q ^ 0.6-1 for the slope of the radiative emissivity. If we adopt q = 0.7 instead of g = 1, our 
predicted 70-160 /im fluxes increase by factors of 2-4. Combining these two modiflcations 
increase our predicted fluxes by a factor of ~ 5-10 at both wavelengths. 

Several observations could check whether these modiflcations of our standard model are 
reasonable. By analogy with the Spitzer Vega data, detection of an outflowing wind of small 
grains in a debris disk around a solar-type star provides a clean measurement of r2 and a 
better constraint on our predicted IR excesses. Measurements of q for larger samples of 
debris disks allows a better assessment of our assumptions for the grain emissivity. 



6. CONCLUSIONS 

Our calculations provide a robust picture for the formation of planets and debris disks 
from a disk of icy planetesimals and set the context for the evolution of dusty debris in a 
dynamic system of planets. The results of this study provide a framework for interpreting 
existing observations of debris disks around 1-3 M0 stars and suggest new observational 
tests of this picture. 

We describe a suite of numerical calculations of planets growing from ensembles of icy 
planetesimals at 30-150 AU in disks around 1-3 Mq stars. Using our hybrid multiannulus 
coagulation code, we solve for the evolution of sizes and orbits of objects with radii of ~ 1 m 
to > 1000 km over the main sequence lifetime of the central star. These results allow us to 
constrain the growth of planets as a function of disk mass, stellar mass, and semimajor axis. 

Debris disk formation is coincident with the formation of a planetary system. All cal- 
culations of icy planet formation at 30-150 AU lead to a coUisional cascade which produces 
copious amounts of dust on timescales of 5-30 Myr. This dust is observable throughout 
the lifetime of the central star. Because we consider a broad range of input parameters, we 
derive the time evolution of (i) dust produced in the collisional cascade and (ii) the IR and 
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submin emission from this dust as a function of disk mass, stellar mass, and time. 

We divide the rest of this section into (i) theoretical considerations, (ii) observable conse- 
quences, and (iii) observational tests. The theoretical considerations build on the highlights 
of icy planet formation in §3.4. Observable consequences of the calculations follow from the 
discussion in 54. The observational tests of the models are described in 55. 



6.1. Theoretical Considerations 

Icy planet formation at 30-150 AU is self-limiting. Starting with a swarm of < 1 km 
planetesimals, runaway growth produces a set of 100-500 km protoplanets. As the 
protoplanets grow, they stir up leftover planetesimals along their orbits. When the 
leftovers reach high e, collisions produce debris instead of mergers. Because proto- 
planets cannot accrete leftovers rapidly, a cascade of destructive collisions grinds the 
leftovers to dust. Poynting- Robertson drag and radiation pressure then remove the 
dust from the disk. 

The maximum sizes of icy planets at 30-150 AU are remarkably independent of initial 
disk mass, stellar mass, and stellar age. For disks with Xm = 1/3 to 3 around 1- 
3 Mq stars with ages t = 0.1 — It^s, the typical planet has r^ax ~ 1750 km and 
iTT-max ~ 0.005 M^. These objects contain < 3%-4% of the initial disk mass. Although 
this result is also independent of the fragmentation parameters, the finite main sequence 
lifetimes of 1-3 M© stars limits the formation of many large planets in the outer disk. 
Thus, the inner disk produces many more Pluto-mass planets than the outer disk 
(Tables [2] and ED. 

For stars close to the main sequence turnoff, stellar lifetimes and the collisional cascade 
limit the mass in solid objects at 30-150 AU. In the inner disk, the collisional cascade 
removes most of the leftover planetesimals before the central star evolves off the main 
sequence. Thus, the typical mass in small objects is ~ 10% of the initial mass at 30- 
40 AU. In the outer disk, smaller collision rates produce a slower cascade. Thus, the 
central star evolves off the main sequence with ~ 50% of the initial mass remaining in 
1-10 km planetesimals at 125-150 AU. 

The collisional cascade produces copious amounts of dust. Dust begins to form during 
the transition from runaway to oligarchic growth [t = 5-10 Myr), peaks when the first 
objects reach their maximum sizes {t = 10-30 Myr), and then slowly declines t > 30- 
50 Myr). The peak mass in 0.001-1 mm (0.001-1 m) particles is ~ 1-2 lunar masses 
(~ 1 M®). Disks with initial masses Xm = 1/3 to 3 reach these peak masses when the 
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age of the star is ~ 10% to 20% of its main sequence lifetime. Because the timescale 
to form dust is short (~ 10-20 Myr), stars are surrounded by large disks of debris at 
30-150 AU throughout their main sequence lifetimes. 

Radiative processes remove large amounts of mass from debris disks. Radiation pres- 
sure produces a radial wind of small particles containing ~ 60% to 70% of the mass 
removed from the disk. Poynting-Robertson drag pulls the rest of the lost mass into 
the inner disk. Because radiation pressure is more important than Poynting-Robertson 
drag when collision rates are large, we expect more wind (inner disk) emission earlier 
(later) in the evolution. 



6.2. Observable Consequences 

We derive clear observational consequences of the collisional cascade. 

» The dusty debris from the collisional cascade is directly observable. For disks around 1- 
3 Mq stars, the maximum fractional dust luminosity of Ld/ L^r ~ 2 x 10~^ is comparable 



to the maximum dust luminosities of known debris disks ( Backman fc Paresce Ill993 



Rieke et al. II2005I : ISu et al. II2006I : iRhee et al. Il2007al ). The dust temperature at the 



inner edge of a 30-150 AU disk scales with the temperature of the central star; thus, 
the predicted 24 /xm excess is very sensitive to the stellar mass. At 70 /im, the predicted 
excesses scale roughly linearly with disk mass and stellar mass. The predicted 160- 
850 yum excesses depend on the disk mass but are nearly independent of the stellar 
mass. 

For systems with little or no emission from terrestrial dust ([5.8]-[8] < 0.1), mid-IR 
color-color diagrams clearly distinguish debris disks around stars of different masses. 
In a [8]-[24] vs [24]-[70] diagram, 2-3 Mq (1-2 Mq) stars have red (blue) [8]-[24] and 
blue (red) [24]-[70] (Fig [20]). In both cases, the color scales with the initial disk mass 
(Fig [21]). Optical colors and spectra generally provide good estimates for stellar mass; 
thus, these diagrams provide good tests of our model predictions. 



6.3. Observational Tests 

We compare our predictions with observations of A-type stars and solar-type stars. 

» For A-type stars, our calculations are the first to explain the observed rise and fall of 
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debris disk fluxes at 24 /iin (Fig. [22l ICurrie et al. ll2008al Jbl). In our picture, the rise in 



debris disk emission corresponds to the transition from runaway growth - when mergers 
of small planetesimals produce larger protoplanets - to oligarchic growth - when the 
collisional cascade begins to grind leftover planetesimals into dust. When oligarchs in 
the inner disk are close to their maximum sizes of ~ 1750 km, the collisional cascade 
produces a maximum in debris disk emission. For a wide range of initial conditions, 
this maximum occurs at 10-20 Myr. As the collisional cascade moves out through 
the disk, smaller collision rates produce less dust which emits at lower temperatures. 
Thus, the 24 /im excess falls with time. The p redicted rate of decline, t ~" with n 



0.6-0 . 8, is close to the ob served rate {n ^ 0.5-1: iGreaves fc Wyatt II2003I : iRieke et al 



20051 : iRhee et al. Il2007af ) 



At longer wavelengths, the maximum excess is larger and lasts longer than at 24 /xm. 
Predicted mid-IR colors also increase slowly with time. Although larger samples of 
A-type stars with 8 /xm photometry would provide a better test of our models, current 
data for the 70 /im excess and the evolution of the [24]-[70] color agree with our predic- 
tions (Fig. [23] - [2^ . For 2-3 M0 stars near the main-sequence turnoff, our calculations 
also yield a clear maxi mum in the 850 /im flu x. Large surveys, such as the proposed 



JCMT Legacy Survey (JMatthews et al. 1120071 ) and submm observations with ALMA 



Herschel and SOFIA, can test this prediction. 

For solar-type stars, our models match observations of most sources. The predicted 
evolution of the 70-160 /im excesses follows the observed rise at 10-100 Myr, the peak 
at ~ 30-100 Myr, and the decline at > 300 Myr. Although ~ 70% (55%) of observed 
debris disks have fluxes that lie within model predictions, our models underpredict 
fluxes for the brightest sources by a factor of 5-10. Fluxes for models with r2 ~ 0.1 /xm 
and g < 0.7 provide better matches to these observations. To guide our choices for r2 
and g, we require spatially resolved images and submm fluxes for these objects. 

For 1-3 M0 stars with ages ~ 0.1-1 Gyr, current data suggest that solar-type stars 
have a larger range of far-IR excesses than A-type stars. In our models, faster debris 
disk evolution around A-type stars produces a smaller dispersion in far-IR excesses and 
colors for stars with ages of 100 Myr to 1 Gyr. Larger samples of debris disks can test 
this prediction in more detail. 

We also consider observations of Vega, the prototypical debris disk. If we adopt mod- 
els with a blowout radius r2 = 10 /im, we can match observations with a standard 
collisional cascade within a broad torus at 80-200 AU. If the torus contains ~ 1-5 M® 
in large objects with r > 1 cm, the cascade can generate (i) the observed ensemble 
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of grains with r ~ 200-300 yum within the torus and (ii) an outflowing wind of small 
grains with r ~ 1-50 fim. 

This conclusion differs from ISu et al. I (l2006l ). who postulate a recent catastrophic 



collision between two large objects as the source of the dusty Vega wind. Although 
the complete destruction of two large icy objects can produce a massive outflowing 
wind, our results suggest that the dusty wind is short-lived and cannot be rapidly 
replenished by the observed population of larger objects. We show that a steady-stat e 
coUisional cascade can explain the Spitzer data (see also iKenyon fc Bromley II2005I ). 
If our interpretation is correct, sensitive observations at 1-10 mm should detect our 
proposed reservoir of larger objects. 



Matching other observations of debris disks requires more realism in our planet formation 
calculations. Adding binary companions and giant planets provides ways to r nodify the evo- 
lution of the coUisional cascade and to i mpose structu re on rings and tori (e.g. IWilner et al. 
2OO2I : iMoro-Martin fc Malhotra I l2005l : lOuillen I I2OO6I ) . Extending the coagulation calcula- 
tions to smaller sizes allows studies of the formation of winds and other large structures. 
Although these calculations have been prohibitively expensive in computing time, rapid ad- 
vances in computing technology will make these additions possible in the next few years. 



Based on the results described here and in iKenyon &: Bromley I (l2004bl ). we conclude 
that debris disks are the inevitable outcome of icy planet formation in a disk of solid objects. 
The basic structures produce d by this model - broad tori and narrow rings of dust that prop- 



Jayawardhana et al. Ill998l : iKalas 



agat e out through the disk (IKenyon fc Bromlev 112004b ) - are consistent with obsery ations 



2005 



Su et al. 



20061 : iFitzgerald et al. Il2007[ l. The 



(e.g. 

model also explains the time evolution of mid-IR colors and fluxes for debris disks around 

A-type and solar-type stars. 
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A. APPENDIX 



A.l. Growth rates 



In standard coa gulation theory, protoplanets accrete material from a swarrn of planetes- 
imals at a rate (e.g., 



Safronov 1 119691 : iLissauer Ill987l : IWetheriU fc Stewart 1119931 ) 



(Al) 



where r is the radius of a planetesimal, VL is the angular frequency of material in the disk, 
V is the random velocity of planetesimals, and v^sc is the escape velocity of the protoplanet. 
The 1 + {vesc/vY term is the gravitational focusing factor. 

To derive the accretion time, we set t = M/M and substitute the orbital period for the 
angular frequency, 

toe (prP/S) [l + {v^Jv)Y' , (A2) 

where p is the mass density of a planetesimal. Throughout runaway growth and the early 
stages of oligarchic growth Vesc/v ^ 1. Because we are interested in the time to produce 
planets with the same r and p in disks with different P and E, we eliminate p and r. Thus, 
the growth time is roughly 

t<x{P/J:){v/Vescy . (A3) 

This equation sets the typical timescale for planet growth in a disk of planetesima ls. If 

E ~ EnXrr^g-^/^ fEa. fl27D1 and v/vesc ~ constant (Fig. 1; Iwetherill fc Stewart I Il993l : 



Goldreich. Lithwick. fc Sarill2004[ ) 



t oc a^x^-^T, 



iv^-i 



(A4) 



This result is close to the t oc o^x^^'^^Eq ^ derived for the formation of the first 1000 km 
object in our calculations (e.g., Eq. (HTl) ). 

To evaluate possible sources for the extra factor of x^'^^ in our derived accretion times, 
we consider the random velocity v of accreted planetesimals. Shorter growth times require 
smaller random velocities. Thus, we consider processes that damp planetesimal velocities. In 
our calculations, collisions and gas drag can reduce w; dynamical friction and viscous stirring 
increase v. At 30-150 AU, gas drag damps rando m velocities ~ 10-20 times more rapidly 
than collisions (JGoldreich. Lithwick. fc Sari 1120041 ). Thus, we ignore collisional damping and 
concentrate on gas drag. 



Rafikov I (120041 ) investigated the dynamics of small planetesimals and growing proto- 
planets in a gaseous nebula. For the early stages of oligarchic growth, the random velocity 
of planetesimals is 



V/Ve 



^ ^gas 



(A5) 
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Substituting this expression into Eq. flA4|) and adopting a constant gas-to-dust ratio, J2gas ~ 
S, we derive 

t ex a^x J^s-i , (A6) 



with 72 = 271 -|- 1. For typical conditions in planetesimal disks, iRafikov I (12004 ) derived 71 ~ 
1/6 to 1/5. Thus, 72 ~ 1.3 to 1.4, close to the exponent of 1.15 derived in our calculations. 

Our treatment of gas drag probably reduces the exponent of Xm in Eq. ( ]A6p from the 
predicted 1.3-1.4 to 1.15. In our simulations, we assume the gas density declines exponen- 
tially on a timescale tgas = 10 Myr. With typical growth times of 20-40 Myr, the gas density 
is ~ 1% to 10% of its initial value when the first 1000 km objects form in the inner disk. 
Thus, gas drag cannot reduce planetesimal random velocities as efficiently as predicted in 
Eq. IA6[ Reducing drag lowers the exponent. With gas depletion timescales ~ 25% to 50% 
of the growth time, we expect an exponent of 72 ~ 1.1 to 1.2, similar to the 72 = 1.15 in our 
calculations. 



A. 2. Radiation from dust 



In the Appendix of iKenyon fc Bromley I (J2004al ). we briefly described our simple algo- 
rithm for the evolution of particles with sizes smaller than the smallest object - r ~ 1 m - 
followed in the multiannulus coagulation code. This algorithm yields the optical depth in 
very small grains ejected from the system and the optical depth in larger grains evolving 
under the influence of collisions and Poynting- Robert son drag. The optical depth in both 
grain populations allows us to derive the time evolution of the disk luminosity and surface 
brightness in bolometric units. Here, we describe the derivation of grain temperature for 
these populations that yields the predicted time evolution of the broadband spectral energy 
distributions of debris disks. 



As in iKenyon fc Bromley I (J2004al ). we divide objects with sizes smaller than ~ 1 m 
into very small grains, small grains, and large grains. In each annulus k of our calculation, 
radiation pressure ejects very small grains with radii between ri and r2. If Pg is the mass 
density of these grains and Mk is the production rate of very small grains in each annulus, 
the very small grains have an integrated optical depth 



r. 



3(A/r2/ri - 1) >^ 
87[pgr2{l - A/rVr^) ~( 




I 



0'b,k+l 



(A7) 



where ab,k is the inner boundary of an annulus centered at a^, hk is the vertical scale height 
in units of the semimajor axis, and vk k is the orbital velocity in annulus k. 
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For small (r = r2 to 1 mm) and large grains (r = 1 mm to 1 m), we derive the optical 
depth Tk in each annulus. To de r ive th e radial surface brightness and total disk luminosity, 



we follow iKenyon &: Hartmann I (119871 ) and derive the amount of stellar radiation absorbed 
by each annulus. We assume a spherical, limb-darkened star with radius i?^, luminosity L^,, 
and limb-darkening parameter eo = 0.6. For a point P at the outer boundary of annulus k 
with height hp above the disk midplane, rays from the star enter the annulus at a scale height 
hin above (below) the midplane. We compute the length / of the path through the disk and 
derive the optical depth along this path as Tp = (//Aafc)rfc, where Aa^ is the width of the 
annulus. The radiation absorbed along this path is e~'^*'/o, where Jq is the flux incident on 
the boundary of the annulus. Numerical integrations over the stellar surface and the vertical 
extent of an annulus yield the amount of flux absorbed by each annulus, which we convert 
to relative surface brightness. A final numerical integration over the radial extent of the disk 
yields the ratio of the disk luminosity to the stellar luminosity, L^/L^. 

To derive the spectral energy distribution of the disk, we make seve ral assumptions. Con 



sistent with observations of scattered light from resolved debris disks (JBackman fc Paresce 



19931 : iLagrange et al. 1120001 ) . we adopt a single albedo u = 0.25 for all grains. For all A, the 
luminosity in scattered light is then uLa/L^; the thermal luminosity emitted by all grains is 
(1 — u)Lii/L^. In each annulus k, we assume grains emit at a temperature Tj /;, where the 
index i refers to discrete bins in grain size. To derive equilibrium temperatures for these 
grains, we assume the grains have an absorption efficiency €„ oc (X/Xq)^ and radiative effi- 
ciency e^ oc (A/Aq)'^. For most grains in our calculations, the grain size is larger than the 
peak wavelength of radiation emitted by the central star. Thus, the grains efficiently absorb 
stellar photons and p = 0. Large grains with r ^ X emit as blackbodies and have g = 0. 
Smaller grains radiate less efficiently and have q ^ 1. 
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Table 1. Grid of Debris Disk Calculations^ 







Stellar Mass 


in Mq 






^m 


1.0*^ 


l.O'^ 


1.5 


2.0 


2.5 


3.0 


0.33 


41 


19 


15 


20 


18 


18 


0.50 


49 


18 


18 


17 


17 


17 


1.00 


49 


15 


22 


17 


15 


15 


2.00 


41 


15 


30 


17 


16 


15 


3.00 


45 


18 


12 


22 


15 


21 




10.00 


10.00 


2.90 


1.22 


0.65 


0.39 



'^Number of independent calculations for each combina- 
tion of M^, Xm 

^32 annulus models at 30-70 AU and at 70-150 AU 

^=64 annulus models at 30-150 AU 



^Main sequence lifetime in Gvr lDemargue et al. I (|2004l ) 
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Table 2. Median number of Plutos at t = tms/^ for disks around 1 M© stars 

30-37 AU 37-45 AU 45-55 AU 55-67 AU 67-82 AU 82-100 AU 100-123 AU 123-146 AU 



0.33 


40 


44 


33 


32 


19 


12 


3 


1 


0.50 


62 


65 


39 


49 


33 


25 


10 


1 


1.00 


111 


110 


73 


73 


47 


55 


26 


5 


2.00 


172 


194 


134 


155 


116 


84 


58 


33 


3.00 


165 


260 


172 


251 


137 


109 


85 


44 
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Table 3. Median number of Plutos at t = tms/3 for disks around 1.5-3 Mq stars 



^m 


30-37 AU 


37-45 AU 


45-55 AU 


55-67 AU 67 


-82 AU 


82-100 AU 


100-123 AU 


123-146 AU 


1.5 Mq 


0.33 


46 


46 


45 


28 


16 


3 


1 


1 


0.50 


69 


85 


68 


48 


53 


25 


2 


1 


1.00 


102 


136 


128 


101 


98 


67 


28 


2 


2.00 


158 


243 


240 


261 


211 


123 


169 


47 


3.00 


201 


239 


301 


295 


381 


165 


198 


97 


2.0 M© 


0.33 


55 


51 


48 


30 


2 


1 


1 


3 


0.50 


115 


87 


85 


51 


20 


10 


1 


2 


1.00 


172 


187 


231 


123 


82 


50 


1 


3 


2.00 


261 


230 


366 


236 


324 


204 


77 


5 


3.00 


259 


262 


398 


264 


446 


295 


171 


48 


2.5 Mq 


0.33 


105 


92 


63 


55 


19 


1 


5 


2 


0.50 


152 


139 


123 


100 


53 


10 


5 


5 


1.00 


164 


174 


234 


198 


125 


90 


7 


6 


2.00 


223 


230 


278 


243 


158 


108 


10 


8 


3.00 


353 


434 


490 


495 


477 


615 


313 


157 


3.0 Mq 


0.33 


133 


121 


62 


36 


6 


5 


1 


1 


0.50 


127 


153 


129 


103 


50 


2 


2 


1 


1.00 


199 


230 


239 


258 


171 


66 


6 


5 


2.00 


224 


353 


376 


392 


342 


265 


173 


5 


3.00 


428 


598 


479 


657 


570 


756 


578 


272 
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Table 4. Predicted Excesses for Disks Around 1 Mq Stars' 



log t (yr) 


log Ld/L^ 


log F24/F24,0 


log 


; Fm/Fro.o 


log Fi6o/-Fl60,0 


log F85o/-F850,0 










^m 


= 0.33 






5.05 


-4.57 


0.000 




0.053 


0.145 


0.090 


5.15 


-4.57 


0.000 




0.053 


0.145 


0.090 


5.25 


-4.63 


0.000 




0.047 


0.131 


0.082 


5.35 


-4.65 


0.000 




0.045 


0.127 


0.080 


5.45 


-4.67 


0.000 




0.043 


0.123 


0.077 








^m 


^ 1.0 






5.05 


-4.18 


0.000 




0.096 


0.255 


0.172 


5.15 


-4.18 


0.000 




0.096 


0.255 


0.172 


5.25 


-4.29 


0.000 




0.090 


0.243 


0.164 


5.35 


-4.31 


0.000 




0.088 


0.238 


0.161 


5.45 


-4.33 


0.000 




0.086 


0.233 


0.158 








^m 


= 3.0 






5.05 


-4.02 


0.000 




0.162 


0.464 


0.357 


5.15 


-4.06 


0.000 




0.150 


0.410 


0.307 


5.25 


-4.10 


0.000 




0.138 


0.355 


0.257 


5.35 


-4.12 


0.000 




0.133 


0.354 


0.248 


5.45 


-4.14 


0.000 




0.127 


0.332 


0.238 



^The electronic version of this paper contains the complete version of this Table. 
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Table 5. Predicted Excesses for Disks Around 1.5 Mp) Stars^ 



log t (yr) 


log Ld/L^ 


log F24/F24,0 


log 


; Fm/Fro.o 


log Fi6o/-Fl60,0 


log F85o/-F850,0 










^m 


= 0.33 






5.05 


-4.40 


0.001 




0.100 


0.202 


0.110 


5.15 


-4.45 


0.001 




0.090 


0.184 


0.099 


5.25 


-4.51 


0.001 




0.078 


0.163 


0.088 


5.35 


-4.52 


0.001 




0.074 


0.157 


0.085 


5.45 


-4.53 


0.001 




0.070 


0.151 


0.082 








^m 


^ 1.0 






5.05 


-4.01 


0.001 




0.162 


0.321 


0.188 


5.15 


-4.19 


0.001 




0.148 


0.299 


0.178 


5.25 


-4.22 


0.001 




0.140 


0.287 


0.170 


5.35 


-4.24 


0.001 




0.133 


0.275 


0.163 


5.45 


-4.24 


0.001 




0.128 


0.266 


0.157 








^m 


= 3.0 






5.05 


-3.73 


0.002 




0.206 


0.406 


0.261 


5.15 


-4.01 


0.002 




0.199 


0.394 


0.251 


5.25 


-4.04 


0.002 




0.192 


0.382 


0.242 


5.35 


-4.05 


0.002 




0.181 


0.361 


0.227 


5.45 


-4.08 


0.002 




0.167 


0.339 


0.211 



^The electronic version of this paper contains the complete version of this Table. 
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Table 6. Predicted Excesses for Disks Around 2.0 Mp. Stars^ 



log t (yr) 


log Ld/L^ 


log F24/F24,0 


log 


; Fm/Fro.o 


log Fi6o/-Fl60,0 


log F85o/-F850,0 










^m 


= 0.33 






5.05 


-4.28 


0.002 




0.128 


0.215 


0.109 


5.15 


-4.34 


0.002 




0.118 


0.201 


0.102 


5.25 


-4.41 


0.002 




0.108 


0.191 


0.096 


5.35 


-4.43 


0.002 




0.103 


0.183 


0.092 


5.45 


-4.46 


0.002 




0.098 


0.175 


0.087 








^m 


^ 1.0 






5.05 


-3.89 


0.003 




0.195 


0.334 


0.185 


5.15 


-4.00 


0.003 




0.186 


0.320 


0.177 


5.25 


4.11 


0.003 




0.177 


0.309 


0.170 


5.35 


-4.16 


0.003 




0.175 


0.307 


0.168 


5.45 


-4.18 


0.003 




0.161 


0.286 


0.156 








^m 


= 3.0 






5.05 


-3.71 


0.005 




0.248 


0.421 


0.250 


5.15 


-3.98 


0.005 




0.239 


0.409 


0.240 


5.25 


-4.02 


0.004 




0.227 


0.389 


0.226 


5.35 


-4.02 


0.004 




0.207 


0.358 


0.205 


5.45 


-4.07 


0.004 




0.192 


0.333 


0.188 



^The electronic version of this paper contains the complete version of this Table. 
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Table 7. Predicted Excesses for Disks Around 2.5 Mp) Stars^ 



log t (yr) 


log Ld/L^ 


log F24/F24,0 


log 


; Fm/Fro.o 


log Fi6o/-Fl60,0 


log F85o/-F850,0 










^m 


= 0.33 






5.05 


-4.18 


0.004 




0.170 


0.259 


0.120 


5.15 


-4.22 


0.004 




0.153 


0.233 


0.108 


5.25 


-4.34 


0.004 




0.137 


0.207 


0.096 


5.35 


-4.36 


0.004 




0.130 


0.197 


0.093 


5.45 


-4.38 


0.004 




0.123 


0.188 


0.086 








^m 


= 1.0 






5.05 


-3.97 


0.008 




0.237 


0.352 


0.181 


5.15 


-4.08 


0.007 




0.219 


0.329 


0.168 


5.25 


-4.09 


0.000 




0.000 


0.000 


0.000 


5.35 


-4.14 


0.006 




0.201 


0.306 


0.154 


5.45 


-4.14 


0.005 




0.187 


0.287 


0.144 








^m 


= 3.0 






5.05 


-3.83 


0.010 




0.285 


0.424 


0.230 


5.15 


-3.95 


0.009 




0.270 


0.403 


0.217 


5.25 


-3.98 


0.009 




0.258 


0.387 


0.206 


5.35 


-4.03 


0.009 




0.248 


0.372 


0.197 


5.45 


-3.48 


0.087 




0.749 


0.887 


0.558 



^The electronic version of this paper contains the complete version of this Table. 
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Table 8. Predicted Excesses for Disks Around 3.0 Mq Stars^ 



log t (yr) 


log Ld/L^ 


log F24/F24,0 


log 


; Fm/Fro.o 


log Fi6o/-Fl60,0 


log F85o/-F850,0 










^m 


= 0.33 






5.05 


-4.08 


0.011 




0.173 


0.229 


0.099 


5.15 


-4.11 


0.010 




0.164 


0.213 


0.091 


5.25 


-4.14 


0.009 




0.152 


0.198 


0.083 


5.35 


-4.16 


0.008 




0.140 


0.183 


0.076 


5.45 


-4.19 


0.007 




0.134 


0.174 


0.072 








^m 


= 1.0 






5.05 


-3.96 


0.015 




0.254 


0.326 


0.150 


5.15 


-4.01 


0.014 




0.239 


0.310 


0.142 


5.25 


-4.04 


0.013 




0.227 


0.297 


0.135 


5.35 


-4.07 


0.012 




0.216 


0.282 


0.127 


5.45 


-4.10 


0.011 




0.204 


0.269 


0.120 








^m 


==3.0 






5.05 


-3.83 


0.020 




0.318 


0.410 


0.199 


5.15 


-3.88 


0.019 




0.301 


0.390 


0.189 


5.25 


-3.92 


0.030 




0.321 


0.389 


0.187 


5.35 


-3.38 


0.113 




0.661 


0.740 


0.416 


5.45 


-3.01 


0.226 




0.934 


0.987 


0.592 



^The electronic version of this paper contains the complete version of this Table. 
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Table 9. Debris disk loci in color-color space 



M, (Mo) 


xo,yo 


x«,y„ 


Xi,y/ 


1.0 


0.00,0.0 


4.0,0.1 


5.00,0.00 


1.5 


1.50,0.0 


4.0,1.5 


4.50,0.25 


2.0 


1.25,0.0 


2.5,2.5 


3.50,0.50 


2.5 


1.00,0.0 


2.0,4.0 


3.00,1.00 


3.0 


1.00,0.0 


1.5,5.0 


2.25,1.75 
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Table 10. Vega Debris Disk ModeP 



Source -^24,* -FVo,* ^leO,* -^850,* F2i^disk FrQ^disk -F'leO.disfe -F850,disfc Afd,l-50 M 

Vega 7.2 0.8 0.16 0.006 1.5 7.0 4.0 0.091 3.0 30'= 



Model r 


7.2 


0.8 


0.16 


0.006 


4.2 


10.0 


3.0 


0.05 


0.9 


0.2 


Model 2"^ 


7.2 


0.8 


0.16 


0.006 


7.2 


23.0 


8.0 


0.11 


3.8 


0.6 


Model 3<= 


7.2 


0.8 


0.16 


0.006 


10.9 


15.0 


4.0 


0.05 


2.4 


0.3 


Model 4^ 


7.2 


0.8 


0.16 


0.006 


15.5 


23.0 


5.0 


0.07 


4.8 


0.8 



'^Fluxes (F) are in units of Jy; dust mass in 1-50 fim particles (M^^i-so) is in units of 10 "^ M^; dust 
production rate (M) is in units of 10^^ g yr~^. 

'^Dust production rate from lSu et al. I ( 2005[ ). Our analysis suggests a smaller dust production rate, 



M >0.3 X 10^1 gyr-i. 

'^Debris disk model with M^ = 2 Mq, a;„ = 1/3, t = 200 Myr 
'^As in note (c) for M^ = 2 Mq, a;,„ = 1, i = 200 Myr 
°As in note (c) for M^ = 2.5 Mq, x™ = 1/2, t = 200 Myr 
fAs in note (c) for M^ = 2.5 Mq, a;™ = 1, i = 200 Myr 
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Fig. 1. — Evolution of a multiannulus coagulation model with S = 

0.18(ai/30 AU)^^/^ g cm~^ at 30-37 AU around a 1 M0 star. Left: median cumu- 
lative mass distribution at times indicated in the legend. Right: median eccentricity 
distributions at t = (light solid line), t = 10 Myr (filled circles), t = 100 Myr (open boxes), 
t = 1 Gyr (filled triangles), and t = 10 Gyr (open diamonds). As large objects grow in the 
disk, they stir up the leftover planetesimals to e ~ 0.1. Disruptive collisions then deplete 
the population of 0.1-10 km planetesimals, which limits the growth of the largest objects. 
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Fig. 2. — Evolution of the radius of the largest object in each annulus for a MMSN disk 
around a 1 Mq star. The number to the right of each set of points indicates log t in yr from 
the start of the calculation. Large objects with r ~ 1000 km form at the inner edge of the 
disk in 10-100 Myr, in the middle of the disk in 0.3-1 Gyr, and at the edge of the disk in 
10 Gyr. 
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Fig. 3. — Median radii of the largest objects at 1 Gyr for disks around a 1 Mq star. The 
numbers to the right of each set of points indicate x^, the disk mass in units of the MMSN. 
Planets form earlier in the inner portions of the most massive disks. Icy planet formation 
saturates at maximum radii r ~ 1500 km. 
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Fig. 4. — Median cumulative mass distributions at 1 Gyr for annuli at 30-37 AU (left panel) 
and at 123-146 AU (right panel) around a 1 M© star. The legend indicates Xm, the initial 
disk mass in units of the scaled MMSN. In the inner disk, many large planets form and the 
collisional cascade removes nearly all of the material in objects with r < 1-10 km. In the 
outer disk, few large planets form; collisions are inefficient at removing material in small 
objects. 
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Fig. 5. — Radius of the largest object in each annulus at 10 Myr and at 100 Myr for a 
MMSN disk around a 1 Mq star. The black points indicate results for calculations with 64 
annuli; the magenta and cyan points show results for calculations with 32 annuli. The good 
agreement between the 32 annulus and 64 annulus calculations shows that planet formation 
is not sensitive to the size of the grid. 
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Fig. 6. — Median production rate of 0.01-1 fim. objects as a function of time for 30-150 AU 
disks around a 1 M© central star. The legend indicates Xm, the initial disk mass in units of 
the scaled MMSN. As large objects grow during the early stages of the evolution, the dust 
production rate declines. Once large objects start to stir up the leftover planetesimals, debris 
production rises rapidly. After dust production peaks, the collision rate and dust production 
decline slowly with time. For all stars, more massive disks eject more material into a wind 
of small particles. 
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Fig. 7. — Median mass in 0.001-1 mm objects as a function of time for 30-150 AU disks 
around a 1 Mq central star. The legend indicates Xm, the initial disk mass in units of the 
scaled MMSN. During runaway growth, the median dust mass is small and roughly constant 
in time. As planet formation propagates through the disk, the dust mass grows with time. 
Once planets form in the outer disk, collisions and Poynting-Robertson drag removes small 
grains from the disk. 
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Fig. 8. — Time evolution of the radius of the largest object at 40 AU and at 100 AU for 
identical disks around 1 Mq (dot-dashed curves), 2 Mq (solid curves), and 3 Mq (triple 
dot-dashed curves) stars. Planets grow faster around more massive stars and in the inner 
disks of all stars. 
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Fig. 9. — Median radius of the largest object at 100 Myr in each annulus for a scaled MMSN 
disk [xm = 1) around 1-3 Mq stars. The number to the right of each set of points indicates 
the stellar mass in M©. At all disk radii, large objects form faster around more massive 
stars. At 30-60 AU, planet formation saturates at radii r ~ 1000-2000 km on relatively 
short timescales, t ~ 100 Myr (see also Eq. (1411) ). At larger disk radii, planets form more 
slowly and do not reach the maximum radius unless the formation time is shorter than the 
main sequence lifetime. 
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Fig. 10. — Median cumulative mass distributions at 100 Myr for planet formation calculations 
at 30-37 AU around 1-3 Mq stars. Left: Results for models with S = 0.18 (ai/30 AU)"^/^ 
g cm^^. The light solid line indicates the initial mass distribution. The dashed (1 Mq, Xm 
= 1), dot-dashed (2 Mq, Xm = 1/2), and heavy solid (3 Mq, Xm = 1/3) lines show median 
results for the same initial conditions. Right Results for models with a scaled surface density 
1), S = 0.18 (ai/30 AU)-3/=^ {MJMq) g cm'^, and different stellar masses (1 Mg 



[Xn 



0- 



dashed line, 2 M©: dot-dashed line, and 3 Mq: heavy solid line). Although more massive 
planets form around more massive stars, the collisional cascade leads to a small dispersion 
in total disk mass at late times. 
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Fig. 11. — Median cumulative mass distributions at t = t^s for annuli at 30-37 AU (left 
panel) and at 123-146 AU (right panel) for identical disks (S = 0.18 (ai/30 AU)"^/^ g cm^^) 
around 1-3 M© stars. The legend indicates the stellar mass in Mq. In the inner disk, many 
large planets form and the collisional cascade removes a large fraction of the material in 
objects with r < 1-10 km. In the outer disk, few large planets form; collisions are inefficient 
at removing material in small objects. 
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Fig. 12. — Median production rate of 0.01-1 /im objects at 30-150 AU as a function of time 
for scaled MMSN (x^ = 1) around 1-3 Mq central stars. The legend indicates the stellar 
mass in Mq. For scaled MMSN, disks around more massive stars eject much more material 
in very small grains at early times (t < 10-100 Myr). In an ensemble of stars with a variety 
of disk masses, there is wide range of dust production rates. 
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Fig. 13. — Median mass in 0.001-1 mm objects as a function of time for scaled MMSN (xm 
= 1) at 30-150 AU around 1-3 Mq central stars. The legend indicates the stellar mass in 
Mq. For t < 1-3 Myr, icy planet formation produces little dust. At 10-100 Myr, the mass 
in small grains is ~ 1 lunar mass for most disks. At late times, the mass in small grains 
slowly declines to currently undetectable levels. 
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Fig. 14. — Time evolution of the median L^/L^^ (dust luminosity relative to the luminosity 
of the central star) for disks surrounding a 1 Mq star. The legend indicates the disk mass 
in units of the MMSN. More massive disks reach larger peak dust luminosities earlier than 
less massive disks. The typical peak dust luminosity is comparable to the dust luminosity 
of the most luminous debris disks associated with solar-type stars. 
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Fig. 15. — As in Fig. [TUfor the median 70 /im excess (left panel) and the median 850 /xm 
excess (right panel). At both wavelengths, dust emission begins to increase at 5-10 Myr. 
Peak dust emission occurs at 30-100 Myr (70 /im) and 100-300 Myr (850 \xm). When the 
central star evolves off the main sequence, the typical excess at 70 /im (850 /xm) is ~ 2-3 
(3-10) times the flux from the central star. 



-90- 



1.5 


1 ' 


n 


1 1 1 1 1 


- 






/ 






-a 1.0 


- 


/ 


\ \ 


- 


o" 




/ 


N \ 




LL. 




/ 


^ \ 150 AU 




r^ 




ft 


\ \ 




IL 




(' 


\ \ 




U) 




/ 


\ \ 




o 






\ \ 




- 0.5 


/ 










jj 




70 AU ^\^ X^^ 


- 




— »_^ 




^v^^ 




0.0 


■ 




1,1, 


" 



^.U 










^- 


1 ' 1 ' 


1.5 








f 


/'' 


""v \ 150 AU 
\ \ 
\ \ 
\ \ 
\ \ 


1.0 






1 






\ \ 
\ \ 
\ \ 
\ \ 
\ \ 
\ \ 
\ 


0.5 




^ 


1 






70 AU \ 

\ 

s 

N 

N 


0.0 


~ 1 








, 


1 , 1 , 1 " 



8 
log Time (yr) 



10 



7 8 

log Time (yr) 



Fig. 16. — Variation of dust excess with disk size. Left panel: time evolution of the median 
70 /im excess for MMSN disks with outer radii of 70 AU (dashed hne) and 150 AU (sohd 
hne). At late times, smaller disks produce smaller IR excesses, it Right panel: as in the left 
panel for the median 850 /im excess. 
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Fig. 17. — Time evolution of the median Ld/ L^, for MMSN disks surrounding 2 Mq stars 
(left panel) and 3 Mq stars (right panel). The legend indicates the initial disk mass in units 
of the scaled MMSN. The typical maximum dust luminosity, Ld/ Li, ~ 10^^ is comparable 
to the dust luminosity of the brightest debris disks around A-type stars. 
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Fig. 18. — Time evolution of median IR excesses for MMSN disks around 1-3 Mq stars. The 
legend in the lower left panel indicates the stellar mass in solar masses for each curve in all 
panels. Lower left panel: 24 /xm excess. Upper left panel: 70 /im excess. Upper right panel: 
160 /im excess. Lower right panel: 850 /im excess. At 24 /im, the peak excess increases 
dramatically with the temperature of the central star. Thus, hotter stars produce much 
larger 24 /im excesses. At longer wavelengths, the magnitude of the excess is correlated with 
the mass of the central star. Roughly independent of stellar mass, the magnitude of the 
excess at 24-850 nm peaks at 10-30 Myr as observed in debris disks around A-type stars 
(jCurrie et al. Il2008al ). At late times, the 160-850 /im excesses for all stars are ~ 3-5 times 
the flux from the stellar photosphere. 
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Fig. 19. — Time evolution of the median [24]-[70] color as a function of time for 1.5 M© 
(dashed line), 2.0 Mq (solid line), and 3.0 (triple dot-dashed line) M© stars. Debris disks 
around lower mass stars have redder [24]-[70] colors than disks around more massive stars. 
For massive stars (> 2 M©), the [24]-[70] color increases slowly throughout the main sequence 
lifetime and then declines just before central star evolves off the main sequence. For lower 
mass stars, the [24]-[70] color reaches a broad maximum at 300 Myr to 1 Gyr and then 
declines. 
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Fig. 20. — Evolution of debris disks with Xm = 1 in color-color space. For icy planet formation 
at 30-150 AU, debris disks around massive stars are hotter than debris disks around less 
massive stars. Thus, debris disks around stars of different masses occupy specific regions of 
the [8]-[24] vs [24]-[70] color-color diagram. 
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Fig. 21. — As in Fig. [20] for disks with different initial masses around 2 Mq stars. Boxes: Xm 
= 1/3, filled circles: Xm = 1, diamonds: Xm = 3. Although more massive disks have redder 
[8]-[24] and [24]-[70], the shape of the color-color track is independent of mass. Thus, the 
color-color diagram isolates stars of different masses. 
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Fig. 22. — Observations of the 24 /xm excess for nearby A-type stars with known ages 



(JRieke et al. I l2005l : ISu et aL I l2006l ) . The hnes show the predicted evolution of the excess 
for debris disk models around 2 M© stars (dot-dashed line: Xm = 1/3; solid line: Xm = 1, 
triple dot-dashed line: x^ = 3). Observations for all but four stars (including one 8 Myr old 
star with log -^24/^24^0 ~ 2) fall within loci defined by our debris disk calculations. Model 
predictions are also consistent with obser vational evidence for a peak in the 24 /xm excess at 
10-20 Myr (see also ICurrie et al. IboOSah . 
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Fig. 23. — As in Fig. [22] for the 70 [iva. excess. Observations for all but 2-3 stars fall within 
the loci defined by the model tracks. Consistent with model predictions, the data suggest a 
peak in the 70 /xm excess at 10-20 Myr. 
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Fig. 24.— Obs ervations of the [24] - [70] color for nearby A-type stars with known ages 
(ISu et al. II2006I ). The hnes show the predicted evolution of the excess for debris disk models 
around 1.5 Mq stars (triple dot-dashed line: Xm = 2) and for 2 Mq stars (dot-dashed line: 
Xm = 1/3; solid line: Xm = 1)- Observations for all but 2-3 stars fall within the model 
predictions. 
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Fig. 25. — Observation s of the 70 /xm excess for nearby solar-type stars with known ages 
(JBeichman et al. II2006I : iHillenbrand et al. II2008I ). The hnes show the predicted evolution 
of the excess for debris disk models around 1 Mq stars (dot-dashed line: x^. = 1/3; solid 
line: x^ = 1, triple dot-dashed line: x^ = 3). Most stars fall within the loci defined by the 
calculations, but many stars are 3-10 times brighter than model predictions. 
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Fig. 26. — As in Fig. [25] for the 160 fim excess. 



